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The  testing,  evaluation,  and  sensitivity  analysis  of  an  in-situ  groundwater  bioreme¬ 
diation  n.odel  is  presented.  BI02D  is  a  two-dimensional,  areal  finite  element  model 
that  incorporates  oxygen  limited  Monod  kinetics.  It  is  validated  using  the  Inter¬ 
national  Groundwater  Modeling  Center’s  Level  1  Testing  Protocol  for  its  ability  to 
model  flow  and  transport  of  a  conservative  pollutant.  Level  2  Testing  is  done  for  four 
hypothetical  cases  involving  the  cleanup  of  a  phenol-contaminated  aquifer.  Predicted 
phenol  concentrations  are  shown  to  be  very  close  and  slightly  greater  than  solutions 
from  the  better  known  and  tested  model  BIOPLUME  II.  It  is  shown  that  BI02D 
outperformed  BIOPLUME  II  in  the  final  test,  which  involved  an  injection-extraction 
well  doublet,  due  to  greater  mass  balance  errors  and  numerical  dispersion  in  BIO- 
PLUME  II.  Two  potential  improvements  to  BI02D  are  presented.  It  is  demonstrated 
that  using  an  iterative  scluti or’  technique  is  preferable  to  the  linearized  method.  The 
proposed  changes  to  the  kinetic  formulation  need  further  testing  however.  A  thor¬ 
ough  sensitivity  analysis  of  BI02D  is  presented.  The  focus  is  on  how  the  predicted 
concentrations  change  with  uncertainty  in  the  often  unknown  biological  parameters. 
It  is  demonstrated  that  individual  sensitivity  alone  is  not  sufficient.  A  technique  of 
combined  sensitivity  called  2*  factorial  design  is  applied  to  these  parameters.  It  is 
shown  that  B102D  is  most  sensitive  to  the  Monod  half  saturation  constant.  A',. 
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Chapter  1 
Introduction 


The  purpose  of  this  research  is  to  test  and  evaluate  an  in-situ  groundwater  bioremedi¬ 
ation  model  named  BI02D.  It  was  selected  as  the  base  model  for  ongoing  optimization 
research  at  Cornell  University.  This  work  will  scrutinize  the  model,  subjecting  it  to 
rigorous  and  thorough  validation  and  comparise  :  testing,  with  the  ultimate  objec¬ 
tives  of  proving  the  model’s  validity,  gaining  valuable  insights  into  its  performance, 
and  observing  its  sensitivity  to  the  often  uncertain  biological  parameters  required  as 
input. 

Motivation  for  Research  It  is  estimated  that  there  are  currently  over  30,000 
contaminated  waste  sites  in  the  U.S.  identified  on  the  Comprehensive  Environmental 
Response,  Compen.'ation  and  Liability  Act  (CERCLA)  list  [Lcgrega  et  al.,  1994]. 
The  Resource  Conservation  and  Recovery  Act  (RCRA)  list  contains  more  than  5000 
additional  hazardous  wMte  facilities.  The  extent  of  this  national  problem  is  enormous 
considering  that  a  significant  number  of  the  7  mi’liot;  undoiground  storage  tanks  in 
the  U.S.  are  estimate-  ,o  be  leaking  [Baker  and  Ilerson,  1994]-  The  U.S.  Army  alone 
has  identified  over  10,000  individual  potentially  contaminated  sites  in  1265  active  and 
closing  installations.  It  is  estimated  that  of  these  over  5000  require  remedial  action 
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[National  Research  Council,  1992]. 

Bioremediation  The  need  to  cleanup  these  sites  has  motivated  the  development 
of  many  new  and  innovative  technologies.  In  1986  CERCLA  was  ammended  to  en¬ 
courage  the  use  of  remediation  technologies  that  would  result  in  permanent  solutions 
[Charbeneau  et  al.,  1992],  One  very  promising  method  that  often  results  in  the  detox¬ 
ification  and  destruction  of  the  contaminant  (rather  than  the  conventional  approach 
of  disposal)  is  bioremediation.  Bioremediation  is  the  use  of  microorganisms  or  mi¬ 
crobial  actions  to  detoxify  and  destroy  environmental  pollutants  [Baker  and  Herson, 
1994]-  Bioremediation  is  not  a  new  tec’  '.ogy;  rather,  it  is  a  new  application  of 
the  same  methods  used  to  treat  and  trans  .  .m  domestic  wastewater  for  the  last  100 
years.  Both  activated  sludge  and  fixed-film  processes  are  based  on  the  exploitation 
of  microorganisms  in  engineered  systems.  The  application  of  bioremediation  to  the 
cleanup  of  hazardous  waste  sites  offers  the  advantages  of: 

•  on-site  remediation,  avoiding  transport  costs  and  liabilities 

•  elimination  of  the  problem  rather  than  merely  moving  it  elsewhere 

•  reduction  of  costs  and  e'eanup  times  in  many  cases 

•  utilizing  an  ecologically  sound  technology  that  applies  natural  proces.ses 

Bioremediation  techniques  have  successfully  been  applied  to  a  wide  range  of  pol¬ 
lutants  including  polyaromatic  hydrocarbons,  volatile  organics  (benzene,  toluene, 
ethylbenzene,  and  xylene  -  often  noted  as  BTEX),  pentachicrophenol,  auid  phenols 
[Charbeneau  et  al.,  1992].  In  addition,  ketones,  esters  and  chlorinated  solvents  (TCE, 
PCE)  have  been  biodegraded  [Baker  and  Herson,  1994]. 

This  work  focuses  on  the  mathematical  modeling  of  in-situ  bioremediation,  or 
biostimulation,  of  polluted  aw^uifers  where  the  existing  subsurface  microorganisms  are 
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stimulated  to  accomplish  the  work  of  degrading  the  pollutants.  This  is  accomplished 
by  controlling  their  environment,  and  providing  an  electron  acceptor  (often  in  the 
form  of  O2  or  )  to  assist  these  organisms  in  ’“eating”  the  waste.  This  application 
of  bioremediation  has  been  successfully  applied  to  numerous  waste  sites  in  the  U.S. 
[U.S.  EPA,  1992a, b  and  Charbeneau  et  ai,  1992].  It  has  been  estimated  that  the  use 
of  in-situ  bioremediation  in  the  cleanup  of  polluted  aquifers  has  reduced  the  cleanup 
time  and  costs  by  a  factor  of  100  in  some  cases  [Lagrtga  et  al.,  199]].  See  Lee  ct  al. 
[1988],  Baker  and  Ilcr.son  [199]],  and  Charbeneau  et  al.  [1992]  for  a  more  complete 
presentation  of  bioremediation. 

1.1  Bioremediation  Modeling 

Groundwater  models  arc  among  the  most  important  scientific  tools  available  for  un¬ 
derstanding  groundwater  processes.  They  have  been  u.sed  extensively  as  planning 
and  management  tools  at  numerous  hazardous  waste  sites  [Sattonal  Rc.teareh  Conn- 
cil,  1992],  Biorcmedialion  models  use  mathematical  representations  of  the  complex 
biodegradation  proccs.ses  and  the  more  widely  understood  flow  and  transport  pro¬ 
cesses  that  govern  groundwater  in  the  subsurface.  Several  bioremediation  models 
have  been  developed  and  presented  in  literature  (see  Section  1.2).  Some  of  these 
attempt  to  represent  the  degradation  in  a  very  thorotigh  and  complete  way,  but  often 
sacrifice  the  ablity  to  apply  the  model  to  field-scale  problems.  Others  employ  greater 
simplifications  of  the  complex  processes  in  order  to  permit  more  efficient  codes  that 
can  be  more  easily  applied  to  real  cleanups. 

BI02D  The  model  studied  in  this  work  is  called  131020.  It  is  a  twodimensional, 
areal  finite  element  model  that  makes  several  assumptions  that  allow  it  to  be  ap¬ 
plied  to  realistic  problems.  It  was  written  by  Dr.  Stewart  Taylor  [Taylor,  1993],  and 
adopted  for  current  optimization  research  at  Cornel!  University.  The  research  cou- 
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pies  differential  dynamic  programming  algorithms  with  BI02D  to  determine  optimal 
treatment  strategies  for  remediation  of  polluted  aquifers.  Sec  Min.sker  and  Shoemaker 
[199^]  And  Culver  and  Shoemaker  [1993]  for  more  details.  BI02D  was  chosen  as  the 
base  simulation  model  chiefly  because  analytical  derivatives  can  be  derived  from  its 
finite  element  formulation. 

BIOPLUME  II  The  most  widely  used  and  recognized  biodegradation  model  is 
BIOPLUME  II  [Rifai  et  al.,  I987bj.  It  is  another  example  of  a  model  that  makes 
several  simplifying  a.ssumptions  in  order  to  obtain  a  ficld-u.seable  code.  BIOPLUME 
II  is  based  on  the  USGS  solute  transport  mode!  .MOC  [h'onikow  and  Bredehoeft, 
1978].  BIOPLUME  II  has  been  successfully  applied  to  at  least  two  field  sites  where 
natural  biodegradation  was  observed,  can  der  Heijde  and  Einawawy,  [1993]  report 
that  BIOPLUME  II  has  received  extensive  verification  and  peer  reviews  of  its  theory. 
An  informal  survey  of  environmental  consultants  that  use  biodegradation  models 
conducted  in  the  summer  of  1993  by  the  author  revealed  that  BIOPLUME  II  is  the 
most  commonly  used  and  recognized.  Although  this  may  be  due  in  part  to  the  paucity 
of  field  scale  biodegradation  models  available,  BIOPLUME  II  is  really  the  industry 
standard  for  biorrmediation  modeling.  BIOPLUME  II  is  used  in  this  work  as  a  basis 
of  comparison  for  bioremediation  testing  of  B102D. 

Model  Validation  Testing  Before  a  groundv/ater  model  can  be  used  as  a  planning 
and  decision  making  tool,  its  cre<Ientials  must  be  established  through  systematic 
testing  and  evaluation  of  the  model’s  performance  (van  der  Heijde  and  Einawawy, 
1992].  The  International  Groundwater  Modtlinf:  Center  (IGWMC)  has  formulated  a 
model  review,  verification  and  validation  procedure der  Heijde  et  ai,  1985]  ihat 
it  recommends  for  all  groundwater  models.  IGWMC  prosents  a  thive-lcvel  testing 
strategy: 
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1.  Level  1:  Verification  by  testing  against  analytical  solutions  for  simple  cases. 

2.  Level  2:  Verification  using  more  complex  cases  or  simplified  real-world  systems 
specifically  designed  to  test  specific  types  of  models  or  code  features. 

3.  Level  3:  Validation  against  independently  obtaine'’  field  or  lab  data. 

Huyakom,  et  al.  [1984]  and  Beljin  [1988]  present  the  applica'  on  of  this  procedure 
to  groundwater  transport  models.  Specific  test  cases  and  specifications  are  recom¬ 
mended  for  Level  1  tests.  Examples  of  Level  2  and  3  tests  are  also  presented.  Beljin 
gives  results  from  model  testing  and  comparisons  of  USGS  Method  of  Characteristics 
(MOC)  model,  Random  Walk,  and  a  finite  element  model  called  SEFTRAN.  This 
work  applies  the  first  two  levels  of  the  IGWMC  protocol  in  the  evaluation  of  BI02D 
drawing  upon  the  work  of  Huyakom,  et  al.  and  Beljin. 

Sensitivity  Analysis  An  important  characteristic  of  a  model  is  its  sensitivity  to 
changes  or  uncertainty  in  the  required  input  parameters.  Biodegradation  models 
typically  involve  the  use  of  several  physical  and  biological  parameters.  This  study 
will  address  the  sensitivity  of  BI02D  to  variations  in  biological  parameters.  The 
overall  objective  will  be  to  discover  which  of  the  parameters  are  most  important  to 
know  with  some  certainty,  and  which  may  be  assumed  from  literature. 

1.2  Literature  Review 

An  excellent  review  of  groundwater  bioremediation  modeling  was  recently  published 
by  Bedicnt  and  Rifai  [1992].  It  contains  a  brief  overview  of  the  bioremediation  tech¬ 
nology  for  contaminated  groundwater,  a  review  of  current  models,  and  an  application 
of  BIOPLUME  II  to  an  aviation  spill  site  in  Michigan.  This  article  is  somewhat  biased 
toward  its  authors’  model,  BIOPLUME  II,  but  it  does  provide  an  excellent  history 
of  the  development  of  biodegradation  models. 
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BI02D  is  an  unpublished  model,  although  Taylor  did  publish  results  of  a  sensi¬ 
tivity  analysis  using  it  in  1993.  Much  of  this  work  is  based  upon  handwritten  notes 
from,  and  personal  conversations  with  Taylor. 

1.3  Specific  Goals  of  this  Research 

The  specific  goals  of  this  work  are: 

1.  To  test  the  ability  of  BI02D  to  correctly  model  the  processes  of  flow  and  trans¬ 
port  of  a  conservative  pollutant.  This  is  presented  in  Chapter  3,  where  IGWMC 
Level  1  testing  is  applied  to  BI02D.  Here  model  predictions  are  compared  to 
analytical  solutions. 

2.  To  test  the  ability  of  BI02D  to  model  the  more  complex  processes  of  biodegra¬ 
dation.  IGWMC  Level  2  testing  is  applied  to  BI02D,  but  here  no  analytical 
solutions  are  available.  Chapter  4  presents  this  work,  where  several  test  cases 
are  evaluated  and  model  predictions  aje  compared  to  BIOPLUME  II. 

3.  To  present  and  evaluate  two  potential  improvements  to  BI02D.  These  are  in¬ 
cluded  as  Chapter  5. 

4.  To  conduct  a  thorough  sensitivity  analysis  of  BI02D  to  the  input  biological 
parameters.  This  includes  the  demonstration  of  a  more  complete  technique 
than  is  commonly  used.  Chapter  6  presents  these  results. 

Chapter  2  will  begin  this  study  with  a  review  of  the  governing  equations  and  a 
presentation  of  the  ways  that  BI02D  and  BIOPLUME  II  solve  them. 


Chapter  2 


Model  Descriptions 

2.1  Introduction 

BI02D  and  BIOPLUME  II  take  very  different  approaches  to  the  problem  of  mathe¬ 
matically  modeling  the  complex  processes  involved  in  biodegradation  of  contaminated 
groundwater.  This  chapter  will  examine  these  two  models  in  some  detail,  looking 
closely  at  how  they  numerically  estimate  the  solution  to  the  governing  partial  dif¬ 
ferential  equations,  and  laying  the  foundation  for  later  comparisons  and  analysis.  A 
summary  of  the  major  differences  between  these  models  is  included. 

BIOPLUME  II  (referred  to  hereafter  as  BIOPLUME)  [Rifai  et  ai,  1987b]  is  a 
relatively  well  known  and  widely  used  model.  The  model  uses  the  finite  difference 
method  to  predict  steady-state  or  transient  hydraulic  heads  and  the  method  of  char¬ 
acteristics  to  predict  contaminant  and  oxygen  concentrations  resulting  from  transport 
and  degradation.  It  is  based  on  the  USGS  solute  transport  model  MOC  [Konikow 
and  Bredehoeft,  1978].  BIOPLUME  models  biodegradation  as  an  instantaneous  re¬ 
action  between  contaminant  and  oxygen.  Microbial  biomass  is  assumed  non-limiting; 
thus,  the  model  does  not  include  biomass  growth  nor  transport. 

BI02D  [Taylor,  1993] h  a  newer  model  that  has  not  been  well  tested  or  validated. 
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It  was  selected  as  a  base  model  for  ongoing  optimization  research  at  Cornell  Univer¬ 
sity.  The  model  predicts  the  steady-state  hydraulic  heads  and  contaminant,  oxygen 
and  microbial  biomass  concentrations  resulting  from  transport  and  degradation.  The 
rate  of  substrate  degradation  is  modeled  using  the  Haldane  variant  of  Monod  kinetics. 

The  two  models  differ  fundamentally  with  respect  to  how  each  solves  the  governing 
equations.  In  addition,  BIOPLUME’s  assumption  of  instantaneous  kinetics  greatly 
simplifies  the  problem  and  in  general  permits  faster  computations. 


2.2  Governing  Equations 

The  two-dimensional,  depth  averaged,  flow  and  transport  equations  for  a  confined 
aquifer  are  as  follows,  where  all  variables  are  defined  in  Appendix  A.  These  are 
the  governing  partial  differential  equations  (PDEs)  that  describe  the  biodegradation 
problem.  Both  BI02D  and  BIOPLUME  approximate  the  solutions  to  them  through 
the  use  of  numerical  techniques. 

Flow  in  a  confined  aquifer: 

=  V-iTVh)  +  Q„  (2.1) 

Transport  of  Substrate: 

=  V  •  (6Z?  VS)  -  bV- vs  -  bn\Rs  +  (Sw  -  S)Qw  (2.2) 
Transport  of  Oxygen: 

Robn^  =  V-{bDVO)  -  bV- VO  -  bnXRo  +  (0,  -  0)Qw  (2.3) 

Transport  of  Biomass: 

d  B 

Ribn—  =  V-{bDVB)  -  bV-VB  ^  bnXRg  +  (B,  -  B)Qw 


(2.4) 
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Darcy’s  Law: 

V  =  -A'V/i  (2.5) 

Hydrodynamic  Dispersion: 

Dijmn  =  +  Do  (2.6) 

In  equations  2.2  -  2.4,  the  left  hand  side  represents  the  time  rate  of  change 
of  substrate,  oxygen  and  biomciss.  The  right  hand  side  represents  the  changes 
due  to  hydrodynamic  dispersion,  convective  transport,  microbial  degradation  and 
fluid  sources/sinks  respectively,  for  substrate,  oxygen  and  biomass.  BI02D  assumes 
steady-state  flow  conditions,  and  the  left  hand  side  of  equation  2.1  equals  zero.  BIO- 
PLUME,  however,  does  not  impose  this  restriction  and  can  thus  solve  transient  flow 
problems. 

Dispersion  Coefficients  In  both  models,  it  is  assumed  that  hydrodynamic  disper¬ 
sion  dom.inates  over  molecular  dispersion,  and  that  the  aquifer  is  isotropic.  Thus  the 
dispersivity  tensors  can  be  defined  in  terms  of  the  longitudinal  and  transverse  dis- 
persivities  of  the  aquifer,  a/,  and  ax.  They  are  related  to  the  dispersion  coefficients 
by 


Dl  =  aiV 

(2.7) 

Dt  —  otV 

(2.8) 

It  follows  that  the  components  of  the  dispersion  coefficent  for  two-dimensional  flow 
in  an  isotropic  aquifer  are  [Bear,  1979]: 

Dzx  — 

(2.9) 

_  l/J 

Dyy  =  aXT?-  +  oiqx 

(2.10) 

Dxy  =  Dyi  =  {ai  -  Qx)-^^. 

(2.11) 
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Table  2.1;  Physical  Parameters  Used  in  BI02D  and  BIOPLUME 


Parameter 

Units 

Description 

n 

(1] 

Effective  porosity.  Assumed  constant  in  BIOPLUME, 

but  n  =  n{i,  y)  in  BI02D. 

Aquifer  thickness,  b  =  b(x,y)  for  both  models. 

R. 

[1] 

Contaminant  retardation  factor;  linear  isotherm: 

fZ,  =  —  consianl 

Kd 

(^1 

Distribution  coefficient. 

T„ 

(¥1 

Transmissivity;  T,*  = 

T»» 

[^] 

Transmissivity;  Tyy  =  A',,  6 

T„  =  T„(x,  y)  ;  Tyy  =  Tyyfx.y)  for  both  models. 

K 

(tI 

Hydraulic  conductivity. 

OL 

[L] 

Longitudinal  dispersivity. 

Qf 

m 

Transverse  dispersivity. 

at  =  at(x,y)  ;  aj  =  arfi.y)  for  both  models. 

S 

(11 

Storativity  of  aquifer.  Required  for  transient  problems 

in  BIOPLUMB  only. 

Retardation  Effects  Both  BI02D  and  BIOPLUME  assume  a  linear  isotherm  for 
substrate  (contaminant)  adsorption.  The  contaminant  retardation  factor  is  calculated 
by 

=  (l  +  •  (2.12) 

Both  models  assume  that  the  oxygen  is  not  retarded.  Microbial  retardation  is  ad¬ 
dressed  in  Section  2.3. 

Physical  Parameters  In  order  to  use  either  BI02D  or  BIOPLUME  several  phys¬ 
ical  parameters  chareicterizing  the  polluted  aquifer  must  be  specified.  These  are 
summarized  in  Table  2.1. 
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Biodegradation  In  equations  2.2  -  2.4,  Xjis,  ^rb  represent  the  rate  of 

change  of  substrate,  oxygen  and  biomass  due  to  biodegradation.  BI02D  and  BIO¬ 
PLUME  differ  in  their  formulation  of  these  degradation  terms.  Since  BIOPLUME 
does  not  model  biomass,  equation  2.4  and  Xrr  are  unique  to  BI02D.  Detailed  ex¬ 
planations  of  these  terms  follows  in  Sections  2.3  and  2.4. 

Boundary  Conditions  Equations  2.1  -  2.4  do  not  completely  define  the  biodegra¬ 
dation  problem;  boundary  and  initial  conditions  are  required  as  well.  The  values 
of  the  initial  concentrations,  boundary  conditions  for  hydraulic  head,  and  appropri¬ 
ate  concentrations  along  the  entire  perimeter  of  the  two-dim  -nsional  domain  must  be 
specified  by  the  user  of  either  model.  The  three  types  of  possible  boundary  conditions 
are: 

1.  Type  1.  Specified  head  or  concentration  boundaries  (Dirichlet  conditions)  for 
which  heads  amd/or  concentrations  are  given. 

2.  Type  2.  Specified  flow  boundaries  (Neumann  conditions)  for  which  the  deriva¬ 
tive  of  the  head  or  concentration  is  given.  A  no-flow  boundary  condition  is  set 
by  specifying  flux  to  be  zero. 

3.  Type  3.  Head  or  concentration  dependent  flow  boundaries  (Cauchy  or  mixed 
boundary  conditions)  for  which  flux  across  the  boundary  is  calculated  given 
a  boundary  head  and/or  concentration  value.  This  type  of  boundary  condi¬ 
tion  is  often  called  mixed  because  it  relates  boundary  heads/concentrations  to 
boundary  flows. 

See  Anderson  and  Woessner  [1992]  or  Bear  [1979]  for  a  more  detailed  description  of 
these  boundary  conditions. 


/ 
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2.3  BI02D 

BI02D  uses  the  Galerkin  finite  element  method  in  space  with  a  variably  weighted 
finite  difference  approximation  in  time  to  find  approximate  solutions  to  equations  2.1 
-  2.4. 

Finite  Element  Method  The  following  is  a  brief  description  of  the  finite  element 
method  (FEM)  used  to  solve  the  problem  of  flow,  transport  and  biodegradation  in 
BI02D.  It  is  not  meant  as  a  complete  coverage  or  derivation  of  the  method.  Rather, 
it  is  meant  to  give  the  reader  a  general  background  of  how  BI02D  uses  FEM  to 
estimate  the  solution  of  the  governing  equations.  Sec  Rtmson  ct  al.  [1971],  Finder 
[1977],  Finder  and  Gray  [1977],  Huyakom  and  Finder  [1983]  and  Segerlind  [1984] 
for  more  complete  coverage  of  FEM  applied  to  groundwater  problems. 

For  convenience  equations  2.1  -  2.4  may  be  redefined  in  linear  differential  operator 
form  as  follows  (see  Appendix  A  for  variable  definitions): 

Li  =  V-(TVh)  +  g„,  (2.13) 

L2  =  R,tn^  -  V-{6£)V5)  +  bV-VS  +  bnXRs  -  (Sw  -  S)Qw  (2.14) 

L3  =  Robn^  -  V-{bDVO)  +  6K-V0  +  bnXjio  -  (0«  -  0)Qw  (2.15) 

Li  =  Rbbn^  -  V-{hDVB)  +  bV-VB  -  bnXjiB  -  (Bw  -  B)Qw  (2.16) 

where  Li  =  L2  L3  =  Li  —  0. 

The  first  step  in  applying  the  FEM  is  to  subdivide  the  domain  into  n/e  elements, 
with  nds  nodes.  It  is  assumed  that  the  solution  to  the  governing  partial  differential 
equations  can  be  approximated  by  the  sum  of  a  series  of  interpolating  basis  func¬ 
tions.  The  objective  is  to  approximate  the  solution  functions  h{x,y,t),  S{x,y,t), 
0{x,y,t),  and  B{x,y,t)  from  equations  2.1  -  2.4.  The  basis  functions  provide  the 
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approximations,  h*,  S*,  0*,  and  B*  as  follows: 


h{x,y,t)  S3  h*(x,y,t)  =  '£,”i\hi(i)wi{x,y)  (2.17) 

S(x,y,t)  S3  S*{x,y,t)  =  '£"i\St(t)wi{x,y)  (2.18) 

0{x,y,t)  S3  0*ix,y,t)  =  Y:”i\0,(t)wi(x,y)  (2.19) 

B{x,y,t)  «  5*(i,y,i)  =  i?,(<)iu,(i,y)  (2.20) 


where  hi{t),  Si{t),  Ot{t),  and  Bi{t)  are  the  values  of  the  hydraulic  head,  and  the 
concentrations  of  substrate,  oxygen  and  biomass  at  node  i,  respectively.  The  term 
Wi{x,y)  is  the  basis  function  for  node  i,  which  is  a  function  of  space  only.  The  terms 
Si{t),  Oi{t),  and  Bi{t)  are  functions  of  time  only. 

The  basis  function  is  an  interpolation  function  equal  to  any  value  between  0 
and  1;  it  will  have  a  value  of  1  if  (x,y)  is  at  node  i,  and  a  value  of  0  for  all  other 
nodes.  wi{x,y)  will  be  between  0  and  1  in  all  elements  that  include  node  i,  and 
0  throughout  any  element  that  does  not  contain  node  i.  There  are  many  possible 
basis  functions  that  meet  this  requirement.  BI02D  uses  linear  basis  functions,  and 
the  nodes  correspond  to  the  corners  of  the  elements.  In  general,  two-dimensional 
elements  may  be  triangular  or  quadrilateral,  but  BI02D  uses  quadrilateral  only. 

Because  h*,  S*,  0*,  and  B*  are  approximations  of  h,  S,  0,  and  B  respectively, 
the  differential  operators  evaluated  using  the  approximate  solutions  define  residual. 


Ri  =  Hh*) 

(2.21) 

R2  =  L2{S*) 

(2.22) 

Ri  =  i3(0*) 

(2.23) 

R4  =  Li{B*) 

(2.24) 

non-zero  errors,  i.e., 
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These  residuals  are  minimized  when  they  are  orthogonal  to  a  set  of  weighting  func¬ 
tions.  In  the  Galerkin  method  the  weighting  functions  arc  chosen  to  be  the  same  as 
the  basis  functions  w,{x,y).  Residual  errors  are  minimized  by  imposing  the  following 
conditions: 


JoW,{x,y)L,(h*)dD  =  0 

(2.25) 

Jow,{x,y)L2{S*)dD  =  0 

(2.26) 

J[jWi(x,y)Li{0*)dD  =  0 

(2.27) 

iom{x,y)U(B*)dD  =  {i 

(2.28) 

for  i  =  1,2,  ...nds  and  where  D  is  the  domain  of  the  region. 

Green’s  theorem  is  applied  to  equalize  the  interelement  continuity  requirement:  for 
the  weighting  and  basis  functions,  and  a  variably  weighted  finite  difference  scheme 
approximates  the  time  derivatives.  The  integrals  may  then  be  evaluated,  and  the 
elemental  coefficients  calculated  element  by  element  and  summed  to  create  the  global 
coefficient  matrices.  The  result  is  the  following  system  of  equations: 

im)  =  +  {n}  (2-29) 

([A/,]  +  Ate  [[iV|  +  E”  1  Qu„,[/’c]''])  {5}'+'  =  (2.30) 

([A/,]  -  At(l  -  e)  [[Af]  +  [Ri]  +  El^i  Qu„,JPcl'])  {sy  + 

E,=1  Qv-,,JPc]’5„  +  +  At(l  -  9){F,y 

([A/o]  +  Ate  [[^]  +  -f  E;^i  =  (2.31) 

([A/„]  -  Ai(l  -  6)  [[yV]  +  [R‘]  +  E™  1  Q^,,[Pc]^)  {Oy  + 

LT=1  QrvJPcYO^  +  Ate{Foy+^  +  At(i  -  e){Foy 
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+  (1  -  6)[Mb\*  +  Me  [[;v]  +  +  Y.T=x  QmAPc]'])  {5}'+'  =  (2.32) 

(0(M4]‘+>  +  (1  -  e)[M,\^  +  Atii  -  0)  [[iV|  +  [/?!]  +  e:1,  Q^JPcY])  {BY  + 

ZZi  QwJPcYB^  +  Ate{Ft,y^^  +  At{\-  emy 

where  {  }  denotes  vectors  and  [  |  denotes  matrices.  The  coefficients  of  these  matrices 
are  defined  in  Appendix  B.  Observe  that  2.30  -  2.32  are  nonlinear  equations  because 
of  [iZ,]  ;  [flo];  and  [A/j]  and  [f?i]. 

Equation  2.29  is  solved  first  for  values  of  the  hydraulic  head  at  each  node.  The 
Darcy  velocities  and  dispersion  coefficients  are  then  calculated  by  2.5  and  2.6.  The 
solutions  for  2.30  -  2.32  are  obtained  by  uncoupling  and  linearizing  these  equations, 
using  a  sequential  procedure.  Equations  2.30  -  2.32  are  solved  in  sequence  with 
matrices  and  vectors  calculated  explicitly  using  the  solution  from  the  previous  time 
step.  Eiquation  2.29  is  resolved  only  when  the  pumping  rates  or  locations  are  changed. 

Taylor  [1993]  reports  that  numerical  testing  of  this  sequential,  linearized  solution 
against  the  fully-coupled,  nonlinear  solution  has  proven  its  accuracy  and  efficiency. 
This  will  be  further  addressed  and  tested  in  Chapter  5. 

Stability  Criteria  The  finite  difference  approximation  in  time  used  by  BI02D  for 
solving  the  solute  transport  equations  is  unconditionally  stable  for  the  implicitness 
factor  (see  Appendix  B),  9  >  0.5.  For  9  <  0.5,  the  solution  is  conditionally  stable. 

It  is  important  to  note  that  in  BI02D  the  user  is  responsible  to  select  the  appropriate 
value  for  9. 

Accuracy  Criteria  The  finite  element  method  is  subject  to  numerical  errors,  in¬ 
cluding  artificial  dispersion  sometimes  called  "numerical  dispersion".  Numerical  er¬ 
rors  often  called  "overshoot”  and  "undershoot"  can  also  occur,  particularly  for  ad- 
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vective  dominated  problems.  Numerical  dispersion  arises  from  errors  associated  with 
the  discretization  of  time  and  space.  To  minimize  these  errors  Anderson  and  Woess- 
ner  [1992]  recommend  that  the  grid  should  be  designed  so  that  the  Peclet  number 
(Pt)  is  less  than  or  equal  to  one.  However,  acceptable  solutions  are  possible  with  F* 
as  high  as  ten  [Huyakom  and  Finder,  1983].  By  aligning  the  element  sides  with  the 
direction  of  flow,  transverse  numerical  dispersion  can  virtually  be  eliminated  [Kinzel- 
barh,  1986].  Imposing  such  Peclet  number  conditions  also  eliminates  ovcrsh(X)t  and 
undershoot.  The  Peclet  condition  is  given  by: 


P*.  =  —  <  1.0 
ft,  =  <  1.0. 

Q, 


(2.33) 

(2.34) 


Similarly  the  time  discretization  should  satisfy  the  accuracy  criteria  that  the 
Courant  number  (Cr)  is  less  than  or  equal  to  one.  The  Courant  condition  is  given 
by: 

Cr  =  ■—  <  1.0.  (2.35) 

Equations  2.33  -  2.35  are  the  accuracy  criteria  for  BI02D. 


Monod  Kinetics  BI02D  models  the  biodegradation  of  substrate  using  a  Haldane 
variant  of  the  Monod  equation  where. 


Ar5  =  BIU- 


Xro  =  BIhF 


5 

K.  +  5  +  ^ 


5 


s 


h,  +  5  +  ^ 
O 


All  variables  are  defined  in  Appendix  A. 


K,  +  0 


0 

(2.36) 

!\»  +  0. 

0 

(2.37) 

.ht  -¥0 

C,Yk, 

(2.38) 
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Microbial  Inhibition  The  inhibition  coefficient  decreases  the  microbial  growth 
rate  at  high  substrate  concentrations,  allowing  toxicity  effects  to  be  modeled  when 
appropriate.  As  A',  — *  oo  or  5  — ♦  0,  inhibition  effects  are  negligible. 

Microbial  Retardation  The  partitioning  of  biomass  between  water  and  solid 
phase  is  critical  to  modeling  the  degradation  of  contaminant  in  the  aquifer.  It  is 
modeled  with  the  Freundlich  isotherm  [Lagrrga  ct  ai,  1994]  given  by: 

B*  =  A'*i#^  (2.39) 

The  corresponding  biomass  retardation  factor  is  defined  by 

/I*  =  +  (2.40) 

wb/re  A't  is  the  Freundlich  constant,  is  the  Freundlich  exponent  and  B*  is  the 
amount  of  bacteria  absorbed  per  unit  weight  of  soil.  =  1  for  a  linear  isotherm. 
Figure  2.1  shows  an  example  of  the  Freundlich  isotherm.  Note  that  iVj  <  1  of*cn  has 
no  physical-chemical  basis. 

Biological  Parameters  The  Monod  kinetics  represented  by  equations  2.36  •  2.38 
include  eleven  biological  parameters  that  must  be  sp>ecified  to  parameterize  B102D. 
These  parameters  can  be  determined  by  collection  and  analysis  of  field  data  or  by 
laboratory  study.  Because  these  parameters  are  site  specific,  and  they  often  vary 
considerably  from  the  lab  to  the  field,  the  values  are  often  not  known  with  much 
certainty.  Chapter  6  will  address  sensitivity  analysis  of  these  parameters.  A  brief 
description  of  each  is  found  at  Table  2.2. 

Assumptions  Several  assumptions  are  made  in  the  formulation  of  BI02D: 

1.  Darcy’s  Law  is  valid  and  hydraulic  head  gra<Jients  are  the  only  driving  mecha¬ 
nism  for  aquifer  flow. 
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Fraundich  l■c(hafm 


Figure  2.1;  Frcundlich  Isotherm  for  Microbial  Retardation 

2.  Steady-state  flow  conditions:  the  model  assumes  that  the  aquifer  heads  respond 
very  quickly  to  pumping  rates,  so  that  the  heads  in  period  T  depend  only  upon 
pumping  within  period  T,  and/or  that  the  time  horizon  for  pumping  is  relatively 
long. 

3.  The  aquifer  is  isotropic  with  respect  to  the  longitudinal  and  transverse  d'sper- 
sivities.  ai  —  QL{x,y)  and  ot  =  Q7-(x,  y)  in  BI02D. 

4.  Ionic  and  molecular  diffusion  are  negligible  as  compared  with  hydrodynamic 
dispersion. 

5.  Oxygen  and  substrate  .are  the  only  rate  limiting  factors  in  microbial  growth; 
nitrogen,  phosphorus  and  other  trace  organics  are  available  in  sufficient  quan¬ 
tities. 
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6.  As  in  Borden  and  Bedient  [1986],  equilibrium  partitioning  of  biomass  between 
the  solid  and  water  phase  is  assumed.  This  effect  is  modeled  in  equations  2.4, 
and  2.36  -  2.38  by  a  retardation  factor,  Ri,,  found  using  a  Freundlich  isotherm 
(equation  2.40). 

7.  There  are  no  changes  in  porosity,  permeability  or  dispersivity  as  a  result  of 
biological  growth  within  the  pore  space.  BI02D  assumes  that  these  physical 
parameters  are  constant  with  time.  See  Molz  et  al.  [1986]  or  Taylor  and  Jaffe 
[1991]  for  more  complex  models  that  consider  these  effects. 

8.  The  complex  subsurfjure  bacterial  population  can  be  represented  as  consisting 
of  a  single  facultative,  heterotrophic  bacterial  type.  In  the  absence  of  hydro¬ 
carbon  they  maintain  a  constant  background  concentration  by  aerobically  or 
anaerobically  consuming  the  background  carbon  found  in  the  aquifer.  In  the 
presence  of  oxygen  and  a  more  easily  degradable  carbon  source  they  preferen¬ 
tially  consume  the  hydrocarbon.  It  is  also  assumed  that  the  yield  coefficient, 
Y  and  death  rate,  r>  are  constant  for  aerobic  and  anaerobic  consumption  of 
substrate  contaminant  auid  background  carbon. 

9.  Background  carbon,  Ce,  remains  constant  in  the  aquifer,  ft  niay  change  form 
through  fermentation,  but  carbon  mass  remains  constant. 

10.  Background  carbon  utilization  is  given  by  a  first  order  rate  constant,  k^. 

11.  Oxygen  is  recharged  at  a  rate  necessary  to  offset  O2  consumed  in  the  aerobic 
degradation  of  background  carbon. 

12.  Oxygen  is  not  retwded. 

13.  Vertical  variations  in  head  and  concentration  of  substrate,  oxygen  and  biomass 
are  negligible,  and  a  depth  averaged  model  is  justified. 
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Table  2.2:  Biological  Parameters  Used  in  BI02D 


ruameter 

Units 

Description 

F 

/2U] 

i 

Stoichiometric  ratio  of  oxygen  to  sabstrate  consumed  in 

degradation.  Determined  from  balaiKed  reaction. 

Neglects  substrate  used  for  net  microbial  growth. 

[3^7) 

Maximum  specific  growth  rate  of  bacteria. 

K. 

m 

Monod  substrate  half-saturation  constant.  Substrate 

concentration  at  one  half  the  maximum  growth  rate. 

Ki 

m 

Substrate  inhibition  constant. 

Accounts  for  substrate  toxicity  effects. 

A'« 

m 

Monod  oxygen  half-saturation  constant.  Oxygen 

concentration  at  one  half  the  maximum  growth  rate. 

Y 

[^) 

Yield  coefBcient.  Defined  as  the  ratio  of  mass 

of  biomess  formed  to  the  mass  of  sabstrate  consumed. 

(si;! 

Endogenous  decay  coeficient.  Accounts  for  energy 

required  for  cell  maintenance,  death  and 

predation  that  reduce  the  growth  rate  of  biomass. 

R* 

(11 

Microbial  retardation  factor,  Freundlich  is  .>tberm: 

^  =  (1  + 

Kk 

(i^^l 

Freundlich  isotherm  constant. 

iV, 

(11 

Isotherm  exponent.  /V»  =  1  for  linear  isotherm. 

G 

m 

Background  carbon  concentration  in  aquifer. 

(e.g.,  humic,  fulvic  acids) 

(jirl 

First  order  decay  rate  of  background  carbon. 
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Table  2.3:  BI02D  Input  Files 


Input  File 

Purpose 

data 

User  specifies  timing  data,  physical  and  biological 

parameters,  well  data,  output  control  and  other  options. 

meah 

User  specifies  spatial  discretisation  data. 

ibc 

User  specifies  boundary  and  initial  conditions  for  problem. 

Input/Output  Files  The  simulation  of  contaminant  transport  and  bioremediation 
in  groundwater  using  BI02D  requires  the  use  of  the  three  input  files  shown  in  Table 
2.3. 
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2,4  BIOPLUME 

BIOPLUME  is  a  two-dimensional,  depth  averaged  biodegradation  model.  It  uses  the 
finite  difference  method  to  predict  hydraulic  heads  and  the  method  of  characteristics 
(MOC)  to  predict  contaminant  and  oxygen  concentrations  resulting  from  transport 
and  degradation.  BIOPLUME  is  built  on  the  MOC  model  [Konikow  and  Bredehoeft, 
1978],  maintaining  the  same  numerical  techniques,  but  .solving  the  transport  equation 
for  both  substrate  and  oxygen. 

Method  of  Characteristics  The  following  is  a  brief  description  of  the  method  of 
characteristics  (MOC)  used  to  solve  the  problem  of  transport  and  biodegradation  in 
BIOPLUME.  It  is  taken  from  Konikow  and  Bredehoefi  [1978]  and  Rifai  et  al.  [1987b]. 
It  is  not  meant  as  a  complete  coverage  or  derivation  of  the  method;  rather,  it  is  meant 
to  give  the  reader  a  general  background  of  how  the  MOC  estimates  the  solution  of 
the  governing  equations.  See  Carder  et  al.  [1964]  Konikow  and  Bredehoefi  [1978] 
for  a  more  complete  description. 

Flow  Equation  BIOPLUME  uses  an  iterative  alternating-direction  implicit  (ADI) 
procedure  to  solve  a  finite  difference  approximation  of  the  flow  equation  (2.1)  [Konikow 
and  Bredehoefi,  1978].  After  the  head  distribution  has  been  completed  for  a  given 
time  step,  the  velocity  of  the  groundwater  flow  is  estimated  at  each  node  using  an 
explicit,  finite  difference  form  of  equation  2.5.  See  Finder  and  Bredehoeft  [1968]  or 
Trescott,  Finder  and  Larson  [1976]  for  a  more  complete  coverage  of  the  ADI  proce¬ 
dure. 

TYansport  Equations  The  technique  used  by  BIOPLUME  is  not  to  solve  partial 
differential  equations  2.2  and  2.3  directly,  but  to  solve  a  nearly  equivalent  system  of 
ordinary  differential  equations  (ODEs).  In  order  to  do  this  the  substantial  or  material 


23 


derivative  is  used.  For  example,  the  substantial  derivative  of  substrate  is  defined  as: 

DS  dS  dS dx  dS dy 
~Di  ~  ~di "" 

The  second  and  third  terms  of  equation  2.41  correspond  to  the  advective  part  of 
equation  2.2.  Thus  equation  2.2  can  be  replaced  by  the  following  system  of  ODEs: 


^  =  14  (2.42) 

^  =  V;  (2.43) 

R,bn^  =  V  {bD  V5)  -  bn\Rs  +  (Sw  -  S)Q„  (2.44) 

The  solutions  of  this  system  of  equations  may  be  given  as; 

X  =  x{i)  ;  y  =  y{t)  ;  S  =  5(t)  (2.45) 

Similarly  equation  2.3  may  be  replaced  by  the  following: 

^  =  14  (2.46) 

S  =  V;  (2.47) 

bn^  =  V-  (bD  VO)  -  bnXfio  +  (Ow  -  0)Qw  (2.48) 

where  solutions  of  this  system  of  equations  may  be  given  as: 

X  =  x{t)  ■,  y  =  y{t) ;  0  =  0(t)  (2.49) 


The  solutions  given  by  equations  2.45  and  2.49  are  called  the  characteristic  curves 
of  partial  differential  equations  2.2  and  2.3  respectively. 

Given  the  solutions  to  equations  2.42  -  2.44  and  2.46  -  2.48,  a  solution  to  the  partial 
differential  equations  (2.2  and  2.3)  may  be  obtained  by  following  the  characteristic 
curves.  This  is  accomplished  numerically  by  introducing  a  set  of  moving  points  or 
particles  that  can  be  traced  within  the  finite  difference  grid.  Every  particle  then 


( 
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corresponds  to  one  characteristic  curve,  and  the  values  of  x,  y,  S  and  0  are  obtained 
as  functions  of  time  for  each  characteristic  [Garder  et  a/.,  1964].  Each  point  has  a 
concentration  and  position  associated  with  it  and  is  moved  within  the  flow  field  in 
proportion  to  the  velocity  at  its  location.  This  may  be  visualized  as  tracing  fluid 
particles  through  a  flow  field  and  noting  changes  in  concentration  as  they  move. 

The  method  of  characteristics  uses  particle  tracking  to  solve  for  changes  in  concen¬ 
tration  due  to  advective  transport,  and  an  explicit  finite  difference  approximation  to 
solve  for  changes  in  concentration  due  to  dispersion,  pumping,  divergence  of  velocity, 
changes  in  saturated  thickness  and  biodegradation. 

The  changes  in  concentrations  due  to  advective  transport  are  calculated  in  BIO¬ 
PLUME  as  follows.  The  first  step  is  to  distribute  a  geometrically  uniform  pattern  of 
traceable  particles  in  the  finite  difference  grid.  The  user  of  BIOPLUME  must  specify 
between  four  and  seven  particles  per  cell,  as  Konikow  and  Bredehoeft  [1978]  found 
this  to  produce  satisfactory  results  in  most  two-dimensional  problems.  In  general, 
using  more  particles  results  in  greater  accuracy,  but  requires  more  computations  and 
longer  run  times.  Each  particle  is  assigned  as  its  initial  concentration  the  initial 
concentration  of  the  node  of  the  cell  cont2uning  the  particle. 

In  e2M:h  time  step  of  the  model  simulation,  all  particles  are  moved  a  distance 
proportional  to  the  size  of  the  time  step  and  the  velocity  at  the  location  of  the 
particle.  The  new  position  of  each  particle  is  computed  using  a  finite  difference 
approximation  of  equations  2.42  -  2.43  and  2.46  -  2.47.  The  x  and  y  velocities  for  any 
particular  particle  are  computed  using  a  bilinear  interpolation  over  the  area  of  half  a 
cell.  After  all  points  have  been  moved,  the  concentration  at  each  node  is  temporarily 
assigned  the  average  of  the  concentrations  of  all  points  then  located  within  the  area 
of  the  cell. 

The  changes  in  concentrations  due  to  hydrodynamic  dispersion,  fluid  sources/sinks, 
divergence  of  the  velocity,  changes  in  the  saturated  thickness  of  the  aquifer  and 
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biodegradation  are  calculated  using  an  explicit  finite  difference  approximation  of 
equation  2.44  and  2.48.  The  resulting  approximate  solution  completes  the  definition 
of  the  characteristic  curves  of  2.2  and  2.3. 

Because  the  processes  of  hydrodynamic  dispersion,  advective  transport,  mixing 
and  biodegradation  are  occuring  continuously  and  simultaneously,  equations  2.42  - 
2.44  and  2.46  -  2.48  should  be  solved  simultaneously.  However,  equations  2.42  -  2.43 
and  2.46  -  2.47  are  solved  by  particle  movement  based  on  implicitly  computed  heads, 
while  equations  2.44  and  2.48  are  solved  explicitly  with  respect  to  concentrations. 
This  source  of  error  is  minimized  using  a  two-step  explicit  procedure  in  which  the 
finite  difference  approximation  of  equations  2.44  and  2.48  are  solved  at  each  node; 
equal  weight  is  given  to  concentration  gradients  computed  at  the  previous  time  level 
and  to  concentration  gradients  computed  from  the  particle  movements. 

Stability  Criteria  The  explicit  numerical  solution  of  the  solute  transport  equa¬ 
tions  have  a  number  of  stability  criteria  associated  with  them  as  discussed  in  Konikow 
and  Bredehoeft  [1978].  These  criteria  restrict  the  allowable  size  of  the  time  step  based 
on  the  grid  size,  dispersion  coefficients,  pumping  rates,  and  calculated  velocities.  In 
BIOPLUME  the  stability  criteria  are  calculated  by  the  program.  The  user  must 
specify  the  overall  time  period  of  interest  and  the  program  determines  the  minimum 
required  time  step  and  number  of  particle  moves  necessary  to  satisfy  stablity. 

Instantaneous  Kinetics  Borden  and  Bedient  [1986]  reported  that  microbial  ki¬ 
netics  had  very  little  effect  on  the  hydrocarbon  (contaminant)  distribution  in  the 
body  of  a  plume,  and  on  the  time  until  contaminant  breakthrough.  They  observed 
that  the  rate  of  oxygen  delivery  is  much  slower  than  the  rate  of  microbial  consump¬ 
tion,  and  that  groundwater  velocity  dominates  the  rate  of  oxygen  availability.  Thus, 
they  replaced  the  original  Monod  kinetics  with  an  instantaneous  reaction  between 
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oxygen  and  substrate.  BIOPLUME  generates  two  sets  of  tracer  particles  (as  de¬ 
scribed  above)  for  tracking  the  change  in  oxygen  and  substrate  with  time.  At  every 
step,  the  solute  transport  equations  (equations  2.2  and  2.3)  are  solved  separately  to 
estimate  the  concentrations  of  oxygen  and  substrate.  The  two  plumes  are  combined 
using  the  principle  of  superposition  to  simulate  the  instantaneous  reaction  between 
them.  In  an  explicit  finite-difference  form,  this  approximation  is  [Rifai  et  al.,  1987b]: 


For  St>^ 

=  St  ~  ^ 

Ot+i  =  0 

(2.50) 

For  St 

St+i  =  0 

0<-n  =  Ot  —  SfF 

(2.51) 

where  S,  0  and  F  are  defined  in  Appendix  A.  Thus,  the  only  biological  parameter 
that  must  be  specified  in  BIOPLUME  is  F. 

BIOPLUME  Options  BI  OPLUME  offers  the  user  three  additional  options  for 
contaminant  transport  and  biodegradation  modeling. 

•  First  order  contaminant  decay. 

•  Biodegradation  from  anaerobic  decay:  modeled  as  first  order  decay  of  contam¬ 
inant. 

•  Biodegradation  from  reaeration;  modeled  as  first  order  decay  of  contaminant. 

In  each  case  the  only  required  input  is  the  coefficient  of  decay.  These  options  are  not 
addressed  nor  tested  in  this  work. 

Assumptions  Konikow  and  Bredehoeft  [1978]  present  the  following  assumptions 
for  the  MOC  m.odel: 
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1.  Darcy’s  Law  is  valid  and  hydraulic  head  gradients  are  the  only  driving  mecha¬ 
nism  for  aquifer  flow. 

2.  The  porosity  and  hydraulic  conductivity  of  the  aquifer  are  constant  with  time, 
and  porosity  is  uniform  in  space. 

3.  Gradients  of  fluid  density,  viscosity  and  temperature  do  not  affect  the  velocity 
distribution. 

4.  Ionic  and  molecular  diffusion  are  negligible  as  compared  with  hydrodynamic 
dispersion. 

5.  Vertical  variations  in  head  and  concentration  of  substrate,  oxygen  and  biomass 
are  negligible. 

6.  The  aquifer  is  homogeneous  and  isotropic  with  respect  to  the  longitudinal  and 
transverse  dispersivities. 

In  addition,  Rifai  ei  ai  [1987bj  \nc\nde  the  following  assumptions  for  BIOPLUME: 

1.  Instantaneous  reaction  between  substrate  and  oxygen  in  the  aquifer,  thus  mi¬ 
crobial  biomass,  and  all  other  nutrients  are  available  in  sufficient  quantities. 

2.  Oxygen  is  not  retarded. 

Required  Parameters  for  MOC/BIOPLUME  In  addition  to  the  physical  pa¬ 
rameters  given  in  Table  2.1,  BIOPLUME  requires  some  additional  parameters  which 
are  associated  with  the  method  of  characteristics  (Table  2.4).  See  Rifai  tt  al.  [1987b] 
and  Konikow  and  Bredehoeft  [1978]  for  details  and  guidance  in  the  use  of  these  vari¬ 
ables. 
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Table  2.4:  Parameters  Required  by  Method  of  Characteristics 


Parameter 

Description 

NPMAX 

Maximum  number  of  particles. 

NITP 

Number  of  iteration  parameters. 

ITMAX 

Maximum  number  of  iterations  in  ADI  procedure 

(used  to  solve  flow  equation). 

TOL 

Convergence  criteria  in  ADI  procedure. 

NPTPND 

Initial  number  of  particles  per  node. 

CELPIS 

Maximum  cell  distance  per  particle  move. 

2.5  Summary  of  Differences  in  Models 

Table  2.5  gives  a  summary  of  the  major  differences  in  how  BI02D  and  BIOPLUME 
approximate  the  solutions  to  the  governing  partial  differential  equations.  An  analysis 
of  how  these  two  models  perform  in  validation  testing  is  the  subject  of  the  following 
two  chapters. 


Table  2.5:  Summary  of  Differences  in  BI02D  and  BIOPLUME 


BI02D 

BIOPLUME 

Procedure  u«ed  to  solve  POEs 

Finite  Element  Method 

Method  of  Characteristics 

Types  of  Flow  Problems  Solved 

steady-state  only 

transient  or  steady-state 

Flexibility  in  Element  Sixes 

great,  can  use  various  sixes  and 

combinations  to  match  problem 

none,  must  use  equally  sixed 

rectangular  cells. 

Kinetics  Used 

Monod 

Instantaneous 

Biolo^csl  Parameters  Req’d 

11  (see  Table  2.2) 

1  (must  specify  F) 

Stability 

unconditionally  stable  for  0  >  .5 

calculated  internally 

Numerical  Errors 

prone  to  numerical  dispersion, 

but  can  be  controlled  by 

finer  discretisation 

minimal  numerical  dispersion 

except  at  wells 

(see  Section  3.7). 

Units  Used 

any  consistent  units 

English 

Additional  Options 

none 

first  order  decay 

degradation  by  reaeralion 

anaerobic  degradation 

Chapter  3 


Validation  Testing 

3.1  Introduction 

The  first  step  in  evaluating  BI02D  was  to  test  its  ability  to  model  solute  transport. 
This  was  done  using  the  International  Ground  Water  Modeling  Center’s  Level  1 
verification  protocol  [Huyakorn  ti  ai,  1984;  Beljin,  1988].  The  objectives  of  this 
process  were  to  check  the  ability  of  its  algorithms  to  solve  the  governing  equations 
of  flow  and  solute  transport  under  a  variety  of  physical  conditions,  and  to  check  the 
accuracy  of  the  computer  code. 

Five  test  problems  with  analytical  or  semi-analytical  solutions  were  used,  ranging 
from  a  simple  one-dimensional  problem  to  a  relatively  complex  problem  involving 
two  pumping  wells.  The  specific  test  problems  applied  in  the  validation  were: 

•  Transport  in  a  semi-infinite  column,  Type  1  inlet  boundary  condition 

•  Transport  from  a  continuous  point  source  in  a  uniform  flow  field 

•  Transport  of  a  .solute  slug  in  a  uniform  flow  field 

•  Transport  in  a  radial  flow  field  caused  by  an  injection  well 
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•  Transport  in  a  non-uniform  flow  field  caused  by  an  injection  -  extraction  well 
doublet 

3.2  Basis  of  Comparison 

In  order  to  evalu'.Le  the  performance  of  BI02D  against  the  analytical  solution,  and 
compare  it  to  BICPLUME,  the  following  criteria  were  used: 

•  Root  Mean  Square  Error,  calculated  as  a  percent  of  the  analytical  solution; 

•  Qualitative  comparisons  such  as  "good”,  "reasonable”,  and  "unacceptable” 
based  on  experience  and  intuition. 


J{observedi  —  analytici)^ 
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3.3  Validation  Test  1:  Transport  in  a 
Semi-Infinite  Column 

The  first  validation  test  considered  one-diniensional  flow  and  transport  in  a  semi¬ 
infinite  column.  A  uniform  flow  field  with  a  constant  concentration  boundary  condi¬ 
tion  was  specified  (Figure  3.1).  The  grid  Peclet  number,  Pt  was  2.0  and  the  initial 
Courant  number,  Cr  was  0.4.  A  complete  problem  statement,  input  specifications, 
the  analytical  solution  and  timing  results  are  included  at  Appendix  C. 


C  = 


Ax  =  10  m.  Ay  =  1  m 


V  = 


m 

dajf 


no  flow 
no  How 

400  m— 


C  =  0 
x 


Figure  3.1:  Schematic  Sketch  of  Validation  Test  1 

Results  from  this  test  are  shown  at  Figures  3.2  -  3.4.  They  indicate  an  acceptable 
match  to  the  analytic  solution  was  attained  using  BI02D  for  both  a  retarded  and 
non-retarded  solute.  The  greatest  source  of  error  was  numerical  dispersion,  which  can 
be  damped  by  reducing  the  time  step  or  by  varying  the  time  weighting  or  implicitness 
factor,  0.  At  the  prescribed  time  step  of  2.5  days,  it  is  necessary  to  vary  the  0  to  0.55 
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to  obtain  a  good  fit  (Figure  3.2). 

It  is  desirable  for  on-going  optimization  research  to  have  the  implicitness  factor 
equal  to  one.  In  this  case  it  is  necessary  to  reduce  the  time  step  to  0.5  day  to 
achieve  comparable  results  (Figure  3.3).  As  shown  in  Table  3.1,  further  refinement 
in  the  solution  can  be  attained  by  reducing  the  time  step  and/or  adjusting  9.  For 
unconditional  stability,  the  implicitness  factor  must  be  greater  than  or  equal  to  0.5. 
Figure  3.4  shows  the  results  for  a  retarded  .solute,  where  the  RMSE  wtis  only  slightly 
greater  than  for  the  non-retarded  case. 


VaidatkjnTest  1:  Case  1  (R  «  1) 


Figure  3.2:  Validation  Test  1,  Case  I  at  25  and  50  days  for  A<  =  2.5  days,  9  =  0.55 
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Validation  T sst  1 :  Casa  1  (R  » 1? 


Figure  3.3:  Validation  Test  1,  Case  1  at  25  and  50  days  for  At  =  0.5  days,  0  =  1.0 


ValdaUon  Tast  1 :  Casa  2  (R  •  2] 


Figure  3.4:  Validation  Test  1,  Case  2  at  25  and  50  days  for  A<  =  2.5  days,  9  =  0.60 
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3.4  Validation  Test  2:  Transport  From  a 
Continuous  Point  Source 

The  second  test  applied  to  BI02D  was  transport  in  a  uniform  flow  field  from  a 
continuous  point  source.  The  continuous  point  source  was  modeled  as  an  injection 
well  with  QCq  =  (Figure  3.5).  The  analytical  solution  assumes  that  the 

source  is  small  enough  that  it  doesn’t  alter  the  flow  field;  this  was  accomolished 
by  specifying  a  small  Q  and  large  Cq.  Fine,  medium  and  coarse  grids  were  used 
to  evaluate  the  model’s  sensitivity  to  spatial  discretization.  A  complete  proble.Ti 
statement,  input  specifications,  the  analytical  solution  and  timing  results  are  included 
at  Appendix  C. 

For  case  1  (Ax  =  60  m,  Ay  =  15  m),  the  best  attainable  match  is  possible  when 
the  implicitness  factor  is  equal  to  0.85  (Figure  3.6),  although  the  solution  with  the 
implicitness  factor  equal  to  1.0  is  only  slightly  worse  (by  0.3%).  At  either  value  BI02D 
produces  acceptable  results  at  1000  and  2000  days.  As  in  Test  1,  numerical  dispersion 
causes  the  greatest  error  and  can  be  damped  by  the  use  of  finer  discretization,  as 
demonstrated  by  case  2  vs.  case  3  (Table  3.2). 

Figure  3.7  shows  the  results  at  2000  days  for  a  transverse  cut  of  the  plume  at  180 
and  900  meters  from  the  source.  The  solution  compares  favorably  to  the  analytical 
solution. 


Distance  (m) 

Figure  3.6:  Validation  Test  2,  Case  1  Centerline  Concentrations  for  Case  1 


I 
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3.5  Validation  Test  3:  Transport  of  a  Solute 
Slug 

This  test  considered  the  transport  of  a  contaminant  slug  released  at  t=0  in  a  uniform 
flow  field.  The  slug  was  modeled  using  an  initial  concentration  at  the  origin  of  3500 
mg  (Figure  3.8).  The  Peclet  number  was  1.25  in  the  x-direction  and  5.00  in  the  y- 
direction.  The  initial  Courant  number  was  0.4.  A  complete  problem  statement,  input 
specifications,  the  analytical  solution  and  timing  results  are  included  at  Appendix  C. 


- -  350  m 

Ax  =  Ay  =  5  m 


X 


Figure  3.8:  Schematic  Sketch  of  Validation  Test  3 

Results  from  Validation  Test  3  show  a  strong  dependence  on  the  time  step  and 
the  implicitness  factor  (Table  3.3).  For  the  specified  time  step  of  one  day,  the  best 
results  are  found  if  0  =  0.65  is  used.  As  shown  in  Figure  3.9  this  solution  contains 
a  13.55%  error.  This  error  is  greatest  in  early  times  and  is  reduced  as  the  plume 
disperses.  If  a  0.1  day  time  step  is  used,  however,  it  is  possible  to  achieve  a  solution 


> 
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with  only  7%  error  and  keep  the  implicitness  factor  at  unity  (Figure  3.10).  As  seen 
in  Figure  3.11,  the  errors  along  two  transverse  sections  of  the  plume  were  much  less 
than  along  the  longitudinal  centerline.  Table  3.3  shows  a  summary  of  results  for  this 
test. 


Validation  Test  3;  Centerline  Concentrations 


Figure  3.9:  Validation  Test  3  Centerline  Concentrations,  A<  =  1  day,  9  =0.65 


43 


VaUdatfon  Tact  3:  S  Oaya  at  x  >  75  m 


VaMatfon  Taat  3;  10  daya  at  x  ••  105  m 


Figure  3.11;  Validation  Test  3  Concentrations  at  x  =  75  m  and  105  m 
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3.6  Validation  Test  4:  Transport  in  a  Radial 
Flow  Field 


The  fourth  validation  test  used  to  evaluate  BI02I)  was  transport  in  a  radially  di¬ 
verging  flow  field  resulting  from  an  injection  well  (Figure  3.12).  The  strength  of  the 
injection  well  was  specified  as  ‘25-?^.  The  resulting  flow  field  is  largest  near  the  well 
and  decreases  with  distance  away  from  the  well.  The  flow  field  and  resulting  trans¬ 
port  generated  by  BI02D  were  checked  against  the  analytical  solution.  The  Peclet 
number  for  this  problem  was  3.33,  and  the  Courant  number  was  0.40.  A  complete 
problem  statement,  input  specifications,  the  analytical  solution  and  timing  results 
are  included  at  Appendix  C. 

Results  from  this  test  demonstrate  the  complexity  of  transport  in  a  non-uniform 
flow  field.  Even  with  a  relatively  fine  spatial  discretization  (Ax  =  Ay  =  1  m)  and  a 
small  time  step  (At  =  1  day)  the  best  results  were  still  16%  in  error  (Figure  3.13). 
In  addition,  the  accuracy  of  the  solution  was  not  significantly  improved  by  reducing 
the  time  step  or  altering  the  implicitness  factor  (Table  3.4).  It  is  interesting  to  note 
that  the  results  from  BIOPLUME  (MOC)  were  not  any  better! 
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<?  =  25^ 


Cross  Section: 


il 
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Validation  Tast  4;  CantarKna  Concentrations 


Figure  3.13:  Validation  Test  4  Radial  Concentrations  at  20  and  40  days 


Table  3.4:  Validation  Test  4  Results,  81020 


IBBESa 

RMSE 

1.0 

1.0 

16.15 

1.0 

0.65 

16.10 

0.5 

1.0 

16.11 

0.1 

1.0 

16.05 
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3.7  Validation  Test  5:  Transport  in  a 

Non-Uniform  Flow  Field  Caused  by  an 
Injection  -  Extraction  Well  Doublet 

The  final  validation  test  used  to  evaluate  BI02D  concerned  solute  transport  between 
a  pair  of  recharging  and  discharging  wells  operating  at  a  constant  flow  rate  (Figure 
3.14).  The  objective  was  to  test  the  ability  of  B102D  to  predict  the  concentration 
breakthrough  curve  at  the  extraction  well.  Both  wells  fully  penetrate  a  uniform 
thickness,  confined  aquifer  that  is  assumed  as  infinite,  homogeneous  and  isotropic. 
The  flow  field  is  assumed  as  steady-state.  Contaminated  water  at  100  ^  is  injected 

3 

at  (60,0)  at  a  flow  rate  of  Q  =  2  and  water  is  pumped  out  at  (150,0)  at  the  same 
rate.  Five  different  cases,  corresponding  to  increasing  dispersion  and  lower  Pe  were 
considered.  In  each,  a  grid  spacing  Ai  =  Ay  =  5  m  was  maintained 

For  the  most  genera!  case  involving  the  combined  influences  of  advection  and  dis¬ 
persion,  an  analytical  solution  does  not  exist  [Huyakom  et  al.,  19S4]-  For  a  more 
limited  case  of  pure  awJvection,  a  semi-analytical  solution  was  devi’ioped  and  pro¬ 
grammed  by  Javandel  et  al.  [1984]-  The  model,  called  RESSQ,  uses  the  complex 
velocity  potential  to  estimate  the  concentration  distribution  in  the  aq"ifer.  It  is  appli¬ 
cable  to  an  aquifer  meeting  the  above  restrictions.  Figure  3.15  shows  the  streamline 
pattern  generated  by  RESSQ.  A  complete  problem  st.atement,  input  specifications, 
and  a  discussion  the  semi-analytical  solution  are  included  at  Appendix  C. 

Figures  3.16  -  3.20  show  the  concentration  breakthrough  curves  at  the  extraction 
well,  as  predicted  by  BI02D  and  BIOPLUME.  The  semi-analytical  solution  obtained 
using  RESSQ  assumes  pure  eidvection.  This  case  is  shown  in  Figure  3.16  where 
ai  =  OT  —  0.01  m  and  the  corresponding  is  500.  Observe  that  the  breakthrough 
curve  predicted  by  BIOPLUME  demonstrated  great  fluctuations  for  which  the  semi- 
analytical  solution  represents  the  approximate  upper  envelope.  These  results  are 
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very  similar  to  those  reported  by  El  Kadi  [1988]  for  MOC.  The  solution  predicted 
by  BI02D  is  quite  poor  as  well.  Numerical  errors  associated  with  such  a  high  Pg  are 
very  visible. 

In  Figures  3.17  -  3.20  where  both  advection  and  dispersion  are  important  observe 
the  following: 

•  The  magnitude  of  the  fluctuations  in  BIOPLUME  become  smaller  as  the  dis¬ 
persion  increases,  and  the  Pg  decreases. 

•  Numerical  errors  in  BI02D  are  no  longer  present  for  Pg  <  5. 

•  It  appears  that  BIOPLUME  demonstrates  more  numerical  dispersion  than 
BI02D,  although  there  is  no  analytical  solution  to  compare  to. 

•  Case  5  (a/,  =  9.1  m,  ap  =  1.8  m)  represents  the  dispersion  in  an  aquifer  that 
will  be  used  extensively  in  Chapter  4.  Observe  in  Figure  3.20  that  BIOPLUME 
appears  to  demonstrate  more  numerical  dispersion  for  this  case. 

These  results  show  the  difficulty  of  accurately  representing  the  transport  of  a 
pollutant  in  a  non-uniform  flow  field  involving  wells.  El  Kadi  considers  this  test 
"severe”,  since  it  involves  mainly  radial  flow  and  curved  flow  lines.  However,  in  the 
application  of  BI02D  or  BIOPLUME  to  an  optimization  model  where  numerous  wells 
and  changing  pumping  rates  are  likely,  this  test  is  an  important  one.  No  RMSE  results 
were  computed  as  the  results  from  RESSQ  did  not  lend  themselves  to  a  meaningful 
comparison. 

The  numerical  problems  of  the  MOC  and  wells  are  documented  elsewhere  [El 
Kadi,  1988]  d  [Konikow  and  Brtdehoeft,  1978],  and  are  only  summarized  here. 
These  errors  ire  a.ssociated  with  the  poor  representation  of  the  radial  flow  near  wells; 
and  by  ‘H"  od  of  estimating  the  solute  mass  removed  from  the  aquifer  at  the 
extraction  iiiJ  introduced  into  the  aquifer  at  the  injection  well.  For  example,  at 
the  extrav,  ..  .  .  /ell  the  MOC  must  regenerate  particles  after  the  particles  representing 
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the  pumped  water  are  removed.  The  net  effect  is  the  artificial  creation  of  mass,  and 
hence,  large  mass  baknce  errors. 

Table  3.5  demonstrates  these  high  mass  balance  errors.  These  errors  were  com¬ 
puted  for  a  1500  day  simulation.  For  Ccise  1  (pure  advection),  the  errors  are  based 
upon  the  RESSQ  semi-anlytical  solution.  For  Cases  2-5  the  errors  were  computed 
based  on  the  total  mass  injected  minus  the  total  mass  extracted.  The  mass  extracted 
was  calculated  by  integrating  the  respective  breakthrough  curve  over  the  1500  day 
simulation.  For  Case  5,  observe  that  BIOPLUME’s  mass  balance  error  was  approxi¬ 
mately  5  times  greater  than  the  error  in  BI02D.  This  will  be  important  to  recall  in 
Chapter  4  when  biodegradation  of  a  similar  case  is  considered. 

The  errors  observed  in  BI02D  for  the  first  two  cases  are  a  function  of  the  very 
high  Pg.  In  cases  involving  greater  dispersion  (and  lower  Pg)  B102D  outperforms 
BIOPLUME  based  on  mass  balance  errors.  In  addition,  it  appears  that  numerical 
dispersion  was  greater  for  BIOPLUME. 

Table  3.5;  Validation  Test  5  Mass  Balance  Errors  (After  1500  Days) 


CASE 

Dispersivities 

Error  in  Percent 

BIOPLUME 

BI02D 

1 

ol  —  —  O  Oi 

25.7 

15.6 

B 

a£  =  ax  =  0.1  m 

18.5 

17.1 

B 

at  =  ax  =  1.0  m 

14.1 

1.7 

B 

at  =  07’  =  5.0  m 

7.9 

0.1 

5 

at  =  9.1  m,  07  =  1.8  m 

5.6 

1.1 

1000 

Time  in  Days 


Figure  3.19:  Validation  Test  5:  07-  =  =  5  m 
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3.8  Summary 

The  application  of  these  recommended  tests  to  BI02D  verify  its  ability  to  satisfac¬ 
torily  solve  the  governing  equations  of  flow  and  solute  transport  under  several  varied 
physical  conditions.  The  accuracy  of  the  computer  code  was  established  as  well.  The 
tests  showed  that  BI02D,  as  is  common  with  Finite  Element  Models  [Anderson  and 
Woessner,  1992],  is  prone  to  numerical  dispersion.  The  user  must  be  careful  to  spec¬ 
ify  a  fine  enough  discretization  to  minimize  this  modeling  error.  This  can  be  assured 
by  following  the  accuracy  criteria  given  by  equations  2.33  -  2.35. 

In  addition,  the  effect  of  At  and  6  are  significant.  In  order  to  keep  9  =  1.0,  as 
is  important  for  some  optimization  applications,  the  time  step  must  be  even  smaller 
than  might  be  necessary  without  this  constraint. 

Table  3.6  shows  how  BI02D  compared  with  BIOPLUME  for  the  first  four  valida¬ 
tion  tests.  BIOPLUME  outperformed  BI02D  by  an  average  of  3%,  largely  because 
the  Method  of  Characteristics  is  not  prone  to  numerical  dispersion  [Anderson  and 
Woessner,  1992;  Garder  et  ai,  /5(?./7except  at  wells. 

Validation  Test  5  demonstrated  the  difficulty  of  accurately  modeling  flow  in  a 
non-uniform  flow  field  involving  well  pairs.  For  this  test  B102D  outperformed  BIO- 
PLUME  based  on  greater  mass  balance  errors  in  the  latter  model. 


Table  3.6:  Validation  Testing  Summary 


RMSE 

Validation  Test 

BIOPLUME 

BI02D 

9 

1 

0.55 

3.08 

0.55 

2 

3.24 

7.68 

0.85 

3 

6.95 

13.55 

0.65 

4 

17.02 

16.10 

0.65 

Average 

6.94 

10.10 

Chapter  4 


Eioremediation  Testing 

4.1  Introduction 

After  validating  BI02D’3  ability  to  accurately  mode!  the  flow  and  transport  of  a 
conservative  contaminant  using  the  IGWMC’s  Level  I  Testing  (Chapter  3),  its  ability 
to  model  the  more  complex  processes  of  bioremediation  was  studied.  This  chapter 
utilizes  the  IGWMC’s  Level  2  Testing  to  accomplish  this.  Unfortunately,  analytical 
solutions  for  oxygen  limited  biodegradation  do  not  yet  exist,  therefore,  the  approach 
taken  was  to  model  four  relatively  simple,  yet  realistic  problems  and  to  compare  the 
solutions  predicted  by  BI02D  and  by  BIOPLUME.  This  resulted  in  some  excellent 
insights  into  these  models.  The  specific  problems  used  were: 

•  Biodegradation  of  a  hydrocarbon  spill 

•  Natural  degradation  of  an  existing  plume 

•  Remediated  cleanup  of  an  existing  plume  with  a  single  injection  well 

•  Remediated  cleanup  of  an  existing  plume  with  an  injection  -  extraction  well 
doublet 
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Biodegradation  Extremes  BI02D  models  biodegradation  as  limited  by  Monod 
kinetics  (Equations  2.12  -  2.14).  In  order  to  evaluate  the  contaminant  concentrations 
predicted  by  BI02D,  it  is  helpful  to  consider  the  bounds  that  it  is  subject  to. 

In  an  instantaneous  model,  degradation  is  limited  only  by  the  amount  of  contami¬ 
nant  and  oxygen  present;  thus,  it  predicts  the  most  optimistic  or  "best  case"  solution 
of  substrate  degradation.  The  concentrations  predicted  by  an  instantaneous  model 
define  a  lower  bound  of  what  one  should  expect  from  a  model  that  limits  degradation. 
BIOPLUME  is  an  instantaneous  model,  and  insomuch  that  it  can  accurately  describe 
flow  and  transport,  defines  a  lower  limit  for  B102D. 

Modeling  of  a  non-degrading  substance  defines  the  most  con.servative  or  pes¬ 
simistic  solution  of  substrate  degradation.  The  concentrations  predicted  for  a  con¬ 
servative  contaminant  define  an  upper  bound  for  BI02D. 

Because  BI02D  predicts  slower  than  instantaneous  kinetics,  but  more  degrada¬ 
tion  than  a  conservative  tracer,  it  should  predict  substrate  concentrations  that  falls 
somewhere  between  these  two  extremes  (Figure  4.1).  As  discussed  in  section  2.3  the 
solution  predicted  by  BI02D  is  a  function  of  the  input  biological  parameters;  the  se¬ 
lection  of  these  parameters  should  change  the  predicted  concentrations,  moving  them 
between  these  upper  and  lower  bounds. 

Test  Aqtufer  The  same  hypothetical  aquifer  was  used  to  test  and  compare  the 
models  for  all  four  degradation  tests.  As  seen  in  Figure  4.2,  the  aquifer  consists  of  a 
two-dimensional,  depth-averaged  confined  ariuifer.  Potential  pumping  locations  are 
indicated  at  nodes  (60,0)  and  (300,0).  A  background  Darcy  velocity  of  0.015  ^  was 
specified;  this  was  the  velocity  reported  by  Borden  and  Bedirnt  [1986]  {ox  the  UCC 
aquifer  in  Conroe.  Texas.  Table  4.1  lists  the  physical  parameters  assumed  for  the 
degradation  tests.  A  spatial  di.scctization  of  =  Aj/  =  5  nj  was  selected  based  on 
the  Accuracy  criteria  giv-n  by  (2..33  -  2.3.5).  (.Several  trial  runs  demonstrated  that  a 
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CcxiMrvaliv*  Solution  (No  Oegradation) 

BI020  Ragion 

BIOPLUME  Solution  (Inslamanaous  Oagradation) 


Figure  4.1:  Biodegradation  Regimes 

Pt^  =  2.78  produced  acceptable  results.)  An  initial  time  step  of  Al  =  1  day  wa.s  se¬ 
lected  based  on  several  trial  runs  and  with  consideration  to  the  order  of  magnitude  of 
the  biological  rate  parameters.  Note  that  for  a  conservative  contaminant,  application 
of  the  Courant  criterion  (2.35)  requires  a  time  step  of  approximately  30  days,  based 
on  the  higher  velocities  found  near  pumping  wells.  Further  analysis  of  the  effects  of 
time  step  size  will  be  addressed. 

The  aquifer  is  axisymetrical  with  respect  to  the  x-axis,  so  it  wa.s  possible  to  model 
the  upper  half  of  the  aquifer  only,  saving  considerably  on  memory  requirements  and 
decreasing  n:n-titne.  Thus,  the  domain  was  discretized  into  1600  rectangular  finite 
elements  and  1701  nodes. 

Boundary  Conditions  The  flow  boundary  conditions  specified  for  the  degradation 
tests  are  Type  1  (constant  head)  along  y  =  0  and  y  =  400  m;  and  Type  2  (no  flow) 
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Table  4.1:  Physical  Parameters  For  Biodegradation  Testing  Problems 


Parameter 

Value 

Darcy  Velocity,  V 

0.015  ^ 

Porosity,  n 

0.29 

Longitudinal  Dispersivity,  ai 

9.1  m 

Transverse  Dispersivity,  oj 

1.8  m 

R» 

1.0 

Aquifer  Thickness,  b 

1.0  m 

Ax 

5.0  m 

Ay 

5.0  m 

At 

1.0  day 

Implicitness  Factor,  9 

1.0 

Peclet  number  (Pcz) 

0.55 

Peclet  number  (Pcy) 

2.78 

Courant  number  (Cr) 

0.003 
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along  X  =  100  m  and  x  =  -100  m  (Figure  4.2).  Boundary  conditions  for  transport 
are: 


C{-60,y,t)  =  0.0 

(4.1) 

0{-60,y,t)  =  3.0 

(4.2) 

B{-60,y,t)  =  0.001  Y 

(4.3) 

Choice  of  Contaminant  The  solution  for  BI02D  is  a  function  not  only  of  the 
code’s  ability  to  model  transport,  but  also  of  the  site-specific  biological  parameters 
(Section  2.3).  For  this  study,  the  hypothetical  pollutant  chosen  Wcis  phenol.  It  is  a 
relatively  common  groundwater  contaminant,  and  a  large  number  of  kinetic  values 
are  reported  in  literature. 

Phenol  is  relatively  mobile  in  soil-water  systems.  It  is  a  common  contaminant 
found  in  leachate  from  hazardous  landfills,  and  is  known  to  be  a  co-carcinogen  [La- 
Grtga  et  ai,  1994].  Based  on  the  intake  of  drinking  water  and  aquatic  organisms, 
the  safe  recommended  level  for  human  health  is  3.5  Several  states,  however,  have 
established  more  stringent  standards  [Oak  Ridge  National  Laboratory,  1989]. 

The  data  required  to  conduct  degradation  testing  and  model  sensitivity  analysis 
were  obtained  from  Taylor  [1993]  who  compiled  parameter  ranges  from  literature,  and 
fit  approximate  probability  distributions  for  each  of  them.  See  Rozich  et  al.  [1983], 
Rozich  et  al.  [1985],  Chang  and  Rittmann  [1987],  Speitel  et  al.  [1987],  Hobson 
and  Mills  [1990],  Lin  [1992],  and  Wagner  [1992]  for  details  of  laboratory  studies 
performed  to  get  these  values.  Specific  biological  input  parameters  used  are  given 
in  Table  4.2.  The  sample  size  refers  to  the  number  of  different  test  results  included, 
several  from  the  same  author.  The  mean  values  were  used  for  biodegradation  testing; 
the  extremes  will  be  used  for  the  sensitivity  analysis  presented  in  Chapter  6.  It  is 
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Table  4.2:  Values  of  Biological  Parameters  Used  in  Biodegradation  Testing 


Sample 

Coefficient  of 

Parameier 

Units 

Size 

Mean 

Maximum 

Mininmm 

Variation 

[^] 

33 

6.48 

15.36 

1.97 

.55 

K. 

lY] 

33 

49.6 

266 

1 

1.32 

Ki 

33 

356.8 

1463 

23 

1.14 

K, 

lY] 

3 

1.0 

2 

0.1 

1 

Y 

(^1 

10 

0.70 

1.02 

0.5 

0.25 

r» 

[3^7! 

3 

0.05 

0.10 

0.001 

1 

A* 

[—’^1 
‘•fnj  J 

3 

15.26 

22.97 

7.55 

0.50 

Nh 

(1) 

3 

1 

1.1 

OJ 

0.10 

WM 

Assumed  F  =  3.0  =  constant 

Assumed  Cc  =  715  =  eonalant 

UPoli 

Assumed  kc  —  10 

■*  =  cmatant 

assumed  that  F,  the  ratio  of  oxygen  to  substrate  used  in  degradation,  is  constant.  It 
is  based  solely  on  stoichiometry,  assuming  complete  oxidation  of  phenol  to  CO2  and 
H2O  from: 


CeHsOH  +  7O2  — ^  6(702  +  3^2^  (4.4) 


Bacterial  Population  It  is  assumed  that  the  complex  subsurface  bacterial  pop¬ 
ulation  can  be  represented  as  a  single  facultative,  heterotrophic  bacteria,  such  as 
Bacillus.  In  the  absence  of  hydrocarbon  they  maintain  a  background  concentration, 
Bi,,  by  aerobically  or  anaerobically  consuming  the  aquifer’s  naturally  occurring  or¬ 
ganic  carbon,  such  as  humic  or  fulvic  acid.  In  the  presence  of  oxygen  and  the  more 
easily  degradable  phenol,  the  microbes  preferentially  consume  phenol.  It  is  assumed 
that  the  yield  coefficient,  Y  and  death  rate,  rj  are  constant  for  aerobic  and  anaerobic 
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consumption  of  phenol  and  background  carbon. 

The  values  for  Cc  and  kc  were  assumed  to  be  constant  at  715  ^  and  10“® 
respectively.  BI02D  assumes  that  the  growth  of  bacteria  on  the  background  carbon 
occurs  at  a  rate  necessary  to  offset  death  and  decay,  so  that 

CckcY  =  (4.5) 

Thus,  for  assumed  values  of  Cc,  k^  Y,  Ri,  and  r^,  the  value  of  Bt,  the  background 
biomass,  is  fixed  by  equation  4.5.  The  corresponding  value  of  must  be  specified  in 
the  model  as  the  initial  condition  for  biomass.  Failure  to  do  this  resulted  in  unstable 
results  in  degradation  test  runs.  From  equation  2.40  (Freundlich  isotherm)  for  the 
assumed  parameters  of  A'j  =  15.26  ^  and  iVj  =  1,  the  microbial  retardation  factor 
Rh  is  100.  Finally,  the  background  biomass  concentration  used  for  this  baseline  case 
from  equation  4.5,  using  the  mean  values  from  Table  4.2,  is  0.001  Assuming  the 
average  cell  weight  of  bacteria  to  be  10~’*^  [Neiderhardt,  1990]^  this  is  equivalent 
to  10®^^.  This  microbial  concentration  fails  in  the  ranges  reported  as  observed  by 
Borden  and  Bedient  [1986]. 

Basis  of  Comparison  In  order  to  compare  the  solutions  generated  by  BI02D  and 
BIOPLUME  the  following  criteria  were  used: 

•  Root  mean  square  error  (when  analytical  solution  is  available;  from  equation 
3.1) 

•  Maximum  substrate  concentration,  Smaz{^),  in  aquifer  at  any  time 

•  Total  substrate  mass  in  aquifer,  5masj(5),  at  any  time 

•  Time  to  achieve  a  1  ppm  standard  {days) 

•  CPU  time  to  complete  simulation  run  (mm) 
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•  Qualitative  comparisons  such  as  "good”,  "reasonable”,  and  "unacceptable” 
based  on  experience  and  intuition. 

In  all  cases,  conservation  of  mass  checks  for  both  substrate  and  oxygen  were  per¬ 
formed. 


65 


4.2  Degradation  Test  1:  Biodegradation  of  a 
Hydrocarbon  Spill 

The  first  degradation  test  considered  the  spillage  of  a  degradable  contami  into 
the  hypothetical  aquifer.  This  is  a  common  scenario  where  leaky  sto;fl^.^r  tanks  dis¬ 
charge  a  hydocarbon  pollutant  into  an  aquifer  for  several  months  or  even  years  before 
discovery.  A  spill  rate  of  QCq  =  5-^  at  (60,0)  (Figure  4.2)  was  assumed.  In  order  to 
make  the  problem  realistic  and  ensure  that  the  assumption  of  an  unchanged  velocity 
field  (as  required  for  the  analytical  solution)  was  met,  a  high  concentration  (16.7 

3 

j')  and  low  flow  rate  (0.00015  ^)  were  specified.  The  biological  parameters  were 
specified  as  the  mean  values  from  Table  4.2. 

Initial  Conditions  The  initial  conditions  are  a  clean  aquifer  (S  =  0  ^),  with 
a  background  oxygen  concentration  of  3.0  ^  and  biomass  concentration  of  0.001 
It  is  assumed  that  the  oxygen  consumed  by  the  facultative  bacteria  in  degrading 
the  background  carbon  is  replenished  from  the  vadose  zone  at  a  rate  necessary  to 
maintain  the  background  concentration. 

Conservative  Tracer  Comparison  For  a  non-degradable  contaminant,  the  ana¬ 
lytical  solution  is  given  by  Equations  C.ll  -  C.15.  Figure  4.3  depicts  that  both  BI02D 
and  BIOPLUME  do  a  good  job  of  predicting  the  solution  for  this  conservative  case. 
The  greatest  error  is  found  in  the  region  closest  to  the  source  because  the  analytical 
solution  at  the  well  is  infinity.  A  time  step  of  100  days  is  sufficient  in  BI02D  to 
achieve  a  good  solution  (5.61%  RMSE),  and  the  improvement  in  reducing  the  time 
step  to  1  day  is  relatively  insignificant  (0.81%  decrease  in  RMSE).  This  small  time 
step  is  necessary  for  the  satisfactory  modeling  of  microbial  degradation,  however. 
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Degradation  Test  1:  Conservative  Tracer 


Figure  4.3:  Degradation  Test  1:  Centerline  Concentrations  for  Conservative  Tracer 
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Degradation  Test  1;  Biodegradation  of  Phenol  Spill 


Distance,  x  (m) 

Figure  4.4:  Degradation  Test  1:  Centerline  Concentrations  for  Phenol  Spill 

Phenol  Degradation  Comparison  Figure  4.4  compares  the  solutions  from  BI02D 
and  BIOPLUME  for  biodegradation  of  the  phenol  spill.  The  most  conservative  so¬ 
lution,  M  depicted  by  the  upper  curve,  represents  no  degradation.  The  BIOPLUME 
solution  represents  a  good  estimate  of  the  instantaneous,  or  most  optimistic  solution. 
The  range  that  one  would  expect  BI02D  solutions  to  fall  is  quite  narrow,  and  the  so¬ 
lution  produced  by  BI02D  (using  mean  values  from  Table  4.2)  is  generally  very  close 
to  BIOPLUME’s  solution.  BI02D  does  overpredict  degradation  near  the  source,  but 
at  the  farthest  extent  of  the  plume  predicts  less  degradation  than  BIOPLUME.  Given 
the  fundamental  differences  in  the  two  models,  and  the  uncertainty  associated  with 
the  biological  input  parameters  for  BI02D,  the  results  are  very  reasonable. 
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4.3  Degradation  Test  2:  Natural  Degradation  of 
an  Existing  Plume 

The  second  degradation  test  considered  the  natural  cleanup  of  the  hydrocarbon  spill 
(modeled  in  Test  1)  in  the  same  hypothetical  aquifer  (Figure  4.2  and  Table  4.1).  It 
is  assumed  that  the  leaky  storage  tanks  were  discovered  and  the  source  of  pollutant 
stopped  after  2000  days.  The  cleanup  is  natural  in  that  no  oxygen  was  injected 
into  the  aquifer:  biodegradation  was  limited  by  the  oxygen  remaining  and  oxygen 
recharged  with  flow  into  the  aquifer.  The  biological  parameters  were  specified  as  the 
mean  values  from  Table  4.2. 

Initial  Conditions  The  oxygen  and  substrate  initial  plumes  were  those  produced 
by  BI02D  for  Test  1  after  2000  days,  and  are  shown  in  Figure  4.5.  The  initial  biomass 
is  assumed  constant  at  a  background  concentration  of  0.001  As  in  Test  1,  it  is 
assumed  that  a  background  concentration  of  oxygen  of  3  ^  is  maintained  everywhere 
(except  in  the  phenol  plume)  by  recharge  with  the  vadose  zone. 

Conservative  Tracer  Comparison  Before  comparing  the  concentrations  pre¬ 
dicted  by  the  two  models  for  natural  degradation,  a  comparison  was  done  of  their 
ability  to  model  the  transport  of  a  identically  shaped  conservative  plume.  Figure  4.6 
shows  how  BI02D  and  BIOPLUME  compare  in  their  abilty  to  model  the  transport 
of  the  plume  given  in  Figure  4.5  if  it  were  a  conservative  tracer.  The  models  pre¬ 
dict  nearly  the  same  tracer  plumes  for  all  of  the  3000  day  simulation.  BIOPLUME 
predicts  a  slightly  higher  maximum  concentration  at  1000  days  (2.8%  higher),  but 
nearly  the  same  at  2000  days  and  identical  at  3000  days.  Observe  that  the  1  ppm 
standard  is  achieved  in  2650  days  without  any  degradation.  This  is  primarily  due 
to  the  hydrodynamic  dispersion  of  the  plume,  and  no  phenol  is  degraded.  Mass  is 
conserved  in  both  models. 


Concentralion  (mgl) 
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□agradation  Tast  2:  Initial  Phenol  Plume 


Figure  4.5:  Degradation  Test  2:  Initial  Conditions 


70 


2.5h 


1  1.5h 


0.5H 


Oagradatlon  Tsst  2:  Corsarvativ*  Pluma 


1000  days 

>  i-. 


BIOPLUME 


BI02D 


/  dr  -  1  day,  thsta  -  1.0 


/  \ 


/ 


«  2000  days 

\ 


;  /  \  3000  days 

/  A  \ 

y  y  \  s 

y  V 


so  100  150  200  250  300  350 

Oistanco.  x  (m) 


400 


Figure  4.6;  Degradation  Test  2;  Centerline  Concentrations  for  Conservative  Tracer 


Phenol  Degradation  Comparison  Figure  4.7  shows  how  the  solutions  for  BI02D 
and  BIOPLUME  compare  for  the  natural  biodegradation  of  the  phenol  plume  given 
in  Figure  4.5.  Comparison  of  Figures  4.6  and  4.7  reveal  the  major  effect  of  natural 
degradation  to  be  at  the  leading  edge  of  the  plume,  where  oxygen  is  not  yet  depleted, 
and  at  the  rear,  where  the  effect  of  recharge  oxygen  is  significant.  The  degraded 
plumes  are  narrower,  but  the  peak  concentrations  are  only  slightly  reduced.  Figure 
4.8  depicts  the  phenol  concentrations  along  transverse  cuts  of  the  plume  after  1000 
days  at  X  =  135  m,  and  after  2000  days  at  x  =  190  m.  Again,  B102D  and  BIO¬ 
PLUME  predict  very  similar  concentration  profiles.  Table  4.3  provides  a  summary 
of  Degradation  Test  2  results.  The  time  to  achieve  a  1  ppm  standard  is  2230  days, 
which  is  420  days  quicker  than  without  any  biodegradation,  at  which  the  total  plume 
mass  remaining  in  the  aquifer  was  less  than  50%  of  the  initial. 
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Table  4.3:  Degradation  Test  2  Results  (Phencl) 


TIME 

MEASURE 

B102D 

BIOPLUME 

Initial  Conditicna 

(^) 

22.92 

22.92 

2508 

2508 

500  days 

Smax 

4.10 

4.23 

2137 

2032 

1000  days 

Sm^s 

2.42 

2.48 

SmA»$ 

1667 

1558 

1500  days 

Smas 

1.63 

1.66 

1269 

1156 

2000  days 

SmAt 

1.16 

1.17 

941 

835 

2500  days 

■Sma^ 

0.84 

0.83 

Sina»» 

675 

576 

3000  days 

Smax 

0.60 

0.59 

Smmtt 

461 

372 

Totals  for 

Time  to  I  ppm  stnd  (days) 

2230 

2230 

3000  Days 

CPU  Time  (mtn) 

146 

4.3 
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4.4  Degradation  Test  3:  Remediated  Cleanup  of 
an  Existing  Plume  With  a  Single  Injection 
Well 

The  third  degradation  test  performed  considered  the  remediated  cleanup  of  the  same 
plume  and  aquifer  used  used  in  Test  2.  It  is  now  assumed  that  the  phenol  spill  was 
discovered  and  a  remediation  strategy  designed.  Specifically  a  fully-penetrating  well 

3 

was  selected  that  injected  oxygenated  water  at  8  ^  at  a  flow  rate  of  4  ^  into  the 
aquifer  at  (60,0)  (see  Figure  4.2).  Again,  the  biological  parameters  were  specified  as 
the  mean  values  from  Table  4.2. 

Initial  Conditions  The  initial  phenol  plume  was  the  same  as  for  Test  2  (Figure 
4.5).  It  is  assumed  that  the  facultative  microbes  maintain  the  assumed  backgound 
concentration  of  0.001  ^  by  aerobic  degradation  of  the  background  carbon,  and  that 
the  O2  used  is  replenished  from  the  vadose  zone. 

Conservative  Tracer  Comparison  The  changes  in  the  velocity  field  caused  by 

.1 

pumping  at  a  rate  of  4  ^  caused  the  plume  to  move  and  disperse  more  quickly 
than  for  Test  2  (no  pumping).  This  is  shown  in  Figure  4.9.  After  1000  days,  the 
plume  has  moved  40  meters  further  downgradient  than  in  Test  2  (Figure  4.6),  and 
the  peak  has  also  been  reduced  from  2.42  to  1.40  Similar  results  are  seen  for 
2000  and  3000  days.  The  time  required  to  meet  the  1  ppm  standard  was  1575  days. 
Thus,  the  effect  of  pumping  is  to  wash  the  phenol  out  and  inc^e^lse  hydrodynamic 
dispersion;  however,  no  phenol  is  degraded.  In  both  models,  there  is  some  error  in 
the  conservation  of  mass  due  to  the  difficulty  of  calculating  the  transport  around  the 
injection  well  [El-Kadi,  I988j.  Table  4.5  shows  that  the  error  decreases  with  time  for 
BI02D,  but  remains  approximately  constant  with  time  in  BIOPLUME. 
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Figure  4.9:  Degradation  Test  3:  Centerline  Concentrations  for  Conservative  Plume 


Table  4.4:  Degradation  Test  3:  Mass  Balance  Errors  for  Conservative  Tracer 


Time  (days) 

BI02D 

Mass  (g)  Error  (%) 

BIOPLUME 

Mass  (g)  Error  (%) 

0 

2508 

0 

2508 

0 

1000 

2613 

4.2 

2539 

1.2 

2000 

2581 

2.9 

2533 

1.0 

3000 

2528 

0.80 

2443 

1.4 

76 


Phenol  Degradation  Comparison  The  significant  effect  of  stimulating  the  biodegra¬ 
dation  of  the  phenol  spill  with  oxygenated  water  is  clearly  seen  in  the  results  of  this 
test.  Table  4.5  provides  a  summary  of  results  for  this  test.  By  injecting  oxygenated 
water  at  a  single  well,  BI02D  predicted  that  the  1  ppm  standard  was  achieved  in  805 
days  (as  compared  to  795  days  predicted  by  BIOPLUME).  This  is  about  one  third 
the  time  that  was  achieved  by  natural  degradation  alone,  and  one  half  that  achieved 
by  pumping  without  oxygen.  Observe  that  of  the  original  2608  grams  of  phenol, 
B102D  predicts  that  only  334  remain  in  the  aquifer  after  1000  days  (as  compared  to 
243  predicted  by  BIOPLUME).  P'or  the  mean  biological  parameters  used,  BI02D  and 
BIOPLUME  produced  very  similar  results.  Figure  4.10  shows  the  centerline  concen¬ 
trations  predicted  by  BI02D  and  BIOPLUME  after  500,  1000  and  1500  days.  Figure 
4.11  shows  the  transverse  concentrations  predicted  by  BI02D  and  BIOPLUME  after 
500  days  at  x  =  140  m  and  after  1000  days  at  x  =  180  m.  The  peak  concentrations 
predicted  are  nearly  identical,  but  BIOPLUME  predicts  more  degradation  at  the 
phenol  plume’s  edges. 

Expanding  the  Time  Step  In  order  to  reduce  the  length  of  the  simulation,  the 
time  step  may  be  increcised.  This  may,  however,  result  in  unacceptable  errors  in 
the  solution.  For  this  degradation  test,  A<  was  expanded  to  2,  3,  4,  5  and  10  days 
to  examine  the  trade  off  between  a  shorter  simulation  and  accuracy.  As  seen  in 
Table  4.6,  the  effect  of  changing  the  time  step  to  2,  3  and  4  days  is  minimal  on  the 
performance  measures,  but  very  significant  in  time  savings.  For  A/  =  2,  3  and  4 
days,  the  predicted  values  for  all  three  measures  of  performance  vary  at  most  by  10 
%,  and  less  than  5%  for  most  measures.  Figure  4.12  shows  the  difference  in  centerline 
concentrations  predicted  by  BI02D  when  the  time  step  is  expanded  from  1  to  4  days. 

In  many  applications,  the  small  loss  in  accuracy  may  be  more  than  compensated  for 
by  the  time  and  money  savings  realized  with  the  use  of  a  larger  At.  Using  a  time 
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Degradation  Test  3:  Phenol  Plume 
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Figure  4.10:  Degradation  Test  3:  Centerline  Concentrations  for  Phenol  Plume 


step  any  larger  than  4  days,  however,  results  in  consistent  errors  in  excess  of  10  %. 

Table  4.7  shows  the  effect  of  varying  the  implicitness  factor,  9.  For  At  =  I  day,  the 
effect  is  negligible.  For  At  =  5  days,  however,  the  solution  is  improved  by  adjusting  9. 
Unfortunately,  no  tested  value  of  9  resulted  in  as  good  of  a  solution  as  was  achieved 
by  the  use  of  a  smaller  time  step. 
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Table  4.5:  Degradation  Test  3  Results  (Phenol) 


TIME 

MEASURE 

BI02D 

BIOPLUME 

Initial  Conditions 

c  (^\ 

22.92 

22.92 

(i?) 

2508 

2508 

500  days 

Smar 

1.74 

1.72 

^mai$ 

879 

756 

1000  days 

^max 

0.70 

0.70 

334 

243 

1500  days 

^max 

0.23 

0.22 

^mast 

84 

41 

Totals  for 

Time  to  1  ppm  stnd  (days) 

805 

795 

1500  days 

CPU  Time  (min) 

72 

40 

Table  4.6:  Time  Step  Expansion  in  Degradation  Test  3 


BI02D  With  At  (days)  of: 

TIME 

MEASURE 

1 

2 

3 

4 

5 

10 

Initial  Conditions 

‘^mar  \  i  ) 

22.92 

22.92 

22.92 

22.92 

22.92 

22.92 

2508 

2508 

2508 

2508 

2508 

2508 

500  days 

Smax 

1.74 

1.76 

1.76 

1.76 

1.64 

1.85 

•Smaif 

879 

895 

903 

907 

807 

988 

1000  days 

Smcx 

0.70 

0.71 

0.72 

0.72 

0.64 

0.80 

Smatt 

334 

345 

351 

353 

295 

416 

1500  days 

Smax 

0.23 

0.24 

0.25 

0.25 

0.19 

0.31 

Sma»f 

84 

89 

93 

94 

67 

128 

Totals  for  Time  to  stnd  (days)  805  812  816  816  765  870 

1500  Days  CPU  Time  (mm)  72  36  24  18  14.5  7.2 


2| 


Degradation  Test  3:  Effect  of  Changing  BioZd  Tlmestep 
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Figure  4.12;  Degradation  Test  3:  Expanding  At 


Table  4.7:  Time  Step  Expansion  and  Various  Values  of  9  in  Degradation  Test  3 


At  =  1  day 

At  =  5  days 

TIME 

MEASURE 

9  =  1.0 

9  =  0.75 

^  =  0.5 

p 

II 

9  =  0.75 

9  =  0.5 

500  days 

Smax 

1.74 

1.75 

1.75 

1.64 

1.81 

1.79 

Smasa 

879 

884 

886 

807 

943 

935 

1000  days 

•S'mar 

0.70 

0.71 

0.71 

0  .64 

0.75 

0.74 

^maaa 

334 

337 

337 

295 

375 

370 

1500  days 

Smax 

0.23 

0.23 

0.23 

0.19 

0.27 

0.26 

Smaaa 

84 

86 

86 

67 

105 

103 

Totals  for 

Time  to  stnd  (days) 

805 

807 

807 

765 

840 

835 

1500  Days 

CPU  Time  (m\n) 

72 

72 

72 

14.5 

14.5 

14.5 
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4.5  Degradation  Test  4:  Remediated  Cleanup  of 
an  Existing  Plume  With  an  Injection  - 
Extraction  Well  Pair 

The  final  degradation  teat  performed  considered  the  remediated  cleanup  of  the  same 
plume  and  aquifer  used  in  Tests  2  and  3.  It  is  now  assumed  that  the  remediation 
strategy  included  a  doublet,  with  one  injection  well,  pumping  oxygenated  water  (8 
into  the  aquifer  at  (60,0),  and  one  extraction  well  pumping  water  out  of  the 
aquifer  at  (300,0).  (See  Figure  4.2).  Both  wells  are  assumed  to  be  fully-penetrating, 

3 

and  pump  at  a  flow  rate  of  2  Thus,  the  total  pumping  effort  remains  the  same 
aa  in  Test  3.  The  biological  parameters  were  specified  eis  the  mean  values  from  Table 
4.2. 

Initial  Conditions  The  initial  phenol  plume  was  the  same  as  for  Tests  2  and  3 
(Figure  4.5).  As  before,  it  is  assumed  that  the  facultative  microbes  maintain  a  con¬ 
stant  background  concentration  of  0.001  ^  by  aerobic  degradation  of  the  background 
carbon,  and  that  the  O2  used  is  recharged  from  the  vadose  zone. 

Conservative  Tracer  Comparison  The  complexities  of  a  non-uniform  flow  field 
caused  by  an  injecection  -  extraction  well  doublet  are  seen  in  the  results  of  this  test. 
Recall  Validation  Test  5,  which  tested  BI02D  and  BIOPLUME  under  a  very  similar 
conditions,  with  a  conservative  contaminant.  In  that  test  we  observed  BIOPLUME’s 
problems  with  wells  that  resulted  in  numerical  dispersion  and  mass  balance  errors. 
Transport  around  the  wells  is  dominated  by  radially  convergent  and  divergent  flow. 

For  a  conservative  tracer  (Figure  4.13)  the  concentration  profiles  predicted  by 
BI02D  and  BIOPLUME  are  somewhat  diff  ’ni.  As  in  Test  3,  the  pumping  results 
in  mass  balance  errors  in  both  models,  although  as  El-Kadi  [1988]  reported,  these 
errors  are  greater  for  the  well  dcmlet  than  for  a  single  well.  As  seen  in  Table  4.8, 
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Degradation  Test  4:  Conservative  Tracer 
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Figure  4.13:  Degradation  Test  4:  Centerline  Concentrations  for  Conservative  Tracer 

BI02D  overpredicted  the  tracer  mass  by  about  6%,  which  remained  approximately 
constant  with  time.  BIOPLUME  overpredicted  the  tracer  mass  by  as  much  as  13.7  %, 
and  this  error  increased  with  time.  In  this  case,  BIOPLUME  actually  demonstrated 
numerical  dispersion.  As  seen  in  Figure  4.13,  the  numerical  dispersion  results  in 
a  lower  peak  concentration  at  500  days.  After  1000  and  1500  days,  however,  the 
additional  mass  artifically  created  by  BIOPLUME  results  in  a  closer  match  of  the 
peaks.  These  errors  are  due  to  BIOPLUME’s  difficulties  in  modeling  transport  near 
wells  [Konikow  and  Bredehoefi,  1978]  and  [El  Kadi,  1988]. 

Phenol  Degradation  Comparison  The  results  for  phenol  degradation  are  seen 
in  Figure  4.15  and  Table  4.9.  Observe  that  BI02D  predicts  greater  degradation, 
or  less  total  mass,  than  BIOPLUME  using  both  Smaz  and  Smass  as  performance 
measures.  In  addition  BI02D  predicts  a  100-day  faster  cleanup.  This  may  be  due  to 


83 


the  fact  that  BIOPLUME  artificially  creates  mass,  as  seen  in  the  conservative  case 
and  in  Validation  Test  5.  In  addition,  the  O2  injected  into  the  aquifer  as  an  electron 
acceptor  is  subject  to  numerical  dispersion,  resulting  in  less  available  at  the  plume  for 
degradation.  This  effect  was  clearly  seen  in  Validation  Test  5  and  is  well  doc^’mented 
as  a  weakness  of  the  MOC  when  dealing  with  wells  [El  Kadi,  1988].  The  combined 
effects  of  these  two  factors  is  more  phenol  mass,  which  is  what  we  observed. 
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Table  4.8:  Degradation  Test  4:  Mass  Balance  Eriors  for  Conservative  Tracer 


Time  (days) 

BI02D 

Mass  (g)  Error  (%) 

BIOPLUME 

Maas  (g)  Error  (%) 

0 

2508 

0 

2508 

0 

500 

2665 

6.2 

2676 

7.0 

1000 

2673 

6.6 

2789 

11.5 

1500 

2675 

6.6 

2819 

13.7 

2.5 
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Degradation  T«st  4:  Phenol  Plume 
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Figure  4.14:  Degradation  Test  4:  Centerline  Concentrations  for  Phenol  Plume 
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Table  4.9:  Degradation  Test  4  Results  (Phenol) 


TIME 

MEASURE 

B102D 

BIOPLUME 

Initial  Conditions 

(^) 

22.92 

22.92 

(ff) 

2508 

2508 

500  day* 

2.15 

2.19 

1138 

1099 

•S^mdx 

0.98 

1.10 

543 

544 

1500  days 

^mar 

0.44 

0.55 

^mafs 

211 

216 

Totals  for 

Time  to  1  ppm  stnd  (days) 

987 

1078 

1500  days 

CPU  Time  (min) 

72 

22 
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4.6  Summary 

In  this  chapter,  four  bioremediation  tests  were  performed  on  BI02D,  with  compar¬ 
isons  made  to  BIOPLUME.  Each  considered  a  different  problem  involving  the  spill 
and  cleanup  of  a  phenol  pollutant  into  a  confined  aquifer.  BI02D  simulations  were 
performed  with  the  eleven  biological  parameters  (Table  4.2)  at  their  mean  values. 

For  the  first  three  tests,  the  two  models  predicted  remarkably  similar  results.  In 
each,  BIOPLUME  predicted  slightly  more  degradation,  which  is  consistent  with  what 
we  expect.  In  Test  4,  BI02D  predicted  slightly  more  degradation  that  BIOPLUME. 
This  is  most  likely  due  to  errors  associated  with  the  MOC  and  wells. 

These  results  can  be  summarized  as  follows: 

•  BI02D  predictions  of  contaminant  concentrations  were  very  reasonable,  falling 
within  the  range  expected  in  all  tests  except  Test  4 

•  Although  both  models  have  difficulty  repre.senting  flow  and  transport  around 
wells,  BI02D  appears  to  have  outperformed  BIOPLUME  in  a  case  involving 
injection  and  extraction 

•  In  general  BIOPLUME  simulations  were  faster  than  BI02D’s.  This  is  due  to  the 
small  time  step  required  by  BI02D  during  the  early  weeks  of  the  cleanup  when 
degradation  is  the  dominant  process.  Chapter  5  will  address  an  improvement 
that  will  permit  much  more  competitive  run  limes. 

•  The  expansion  of  the  Af  from  1  day  to  2-4  days  was  found  to  be  quite  reasonable 
for  Test  3.  The  loss  of  accuracy  was  relatively  small  and  the  run  times  were 
considerably  faster.  For  many  applications  a  larger  time  step  may  be  justified. 

•  The  adjustment  of  0,  the  time  implicitness  factor  was  found  to  be  insignificant 
to  solution  accuracy.  This  is  due  to  the  relatively  small  time  step  required  by 
the  degradation.  For  most  applications  taking  B  =  1.0  is  satisfactory. 


Chapter  5 


BI02D  Program  Improvements 

The  previous  chapters  have  presented  and  evaluated  the  biodegradation  model  BI02D 
as  written  by  Taylor.  It  is  the  purpose  of  this  chapter  to  propose  and  test  two 
modifications  lo  the  model  that  may  improve  performance  and  better  represent  the 
processes  involved  in  groundwater  biodegradation. 

5.1  Iterative  Procedure 

As  discussed  in  Section  2.3,  the  solution  technique  for  finding  the  approximate  solu¬ 
tions  to  the  governing  partial  differential  equations  (2.2  -  2.4)  was  to  uncouple  and 
linearize  them.  This  procedure  has  the  advantage  of  being  computationally  simple, 
but  its  accuracy  must  be  further  evaluated.  An  alternative  technique  is  to  solve  these 
equations  by  the  use  of  an  iterative  method.  This  method  is  presented  as  a  way  of 
checking  the  accuracy  of  a  linearized  approach,  and  as  a  potential  refinement  in  the 
solution  procedure. 
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Method  of  Solution  Considering  only  substrate,  equation  2.30  could  be  reformu¬ 
lated  as: 

(M  +  [A|'+')  {<;}'+'  =  ([5,1  +  |fl.l')  {■?}'  + li),l  (5.1) 

where  [i?,]*"*"*  is  a  function  of  {5}^'*"'.  In  the  linearized  method,  is  approx¬ 

imated  by  [/?,]*  and  the  equation  is  solved  in  one  step.  In  an  iterative  method, 
however,  is  updated  after  is  found,  and  the  equation  is  solved  again. 

This  procedure  is  continued  until  the  solution  converges  to  some  acceptable  tolerence. 
This  can  be  represented  by: 

(M,i  +  M!:l){5r' =  (|5.l  + W){s}'  +  li),|  (5.2) 

where  i  =  the  iteration  number.  This  method  is  applied  to  the  partial  differential 
equations  for  substrate,  oxygen  and  biomass. 

The  euclidean  length,  I2  norm,  was  used  as  the  basis  for  convergence  of  the  iter¬ 
ative  method.  It  is  defined  by: 

||x||  =  /ixiF+lxap +  ...  + li„P  (5.3) 

The  residual  errors  in  substrate,  oxygen  and  biomass  convergence  are  given  by: 

-  IIS-!II 

- i5i!!ir~  - 

iis;*i-b;+'ii 

The  iterations  continue  until  the  following  criteria  are  met: 

6s  <0  (5.7) 

60  <0  (5.8) 

6b  <0  (5.9) 


Ml 

'Xj 

I 
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where  l3  is  the  convergence  tolerence.  As  a  matter  of  practicality,  a  maximum  number 
of  iterations,  mxiir  can  be  set  for  the  procedure. 

The  application  of  this  procedure  was  made  to  Degradation  Test  3,  where  a  pol¬ 
luted  aquifer  is  cleaned  up  using  a  single  Injection  well  (see  Section  4.4).  For  this 
test,  the  solutions  predicted  by  BI02D  and  BIOPLUME  were  nearly  identical. 

For  the  baseline  case  (biological  parameters  assumed  as  the  mean  values  from 
Table  4.2),  a  plot  of  residuals  (equations  5.4  -  5.6)  versus  time  is  given  by  Figure  5.1. 
These  residuals  represent  the  error  in  the  linearized  solution.  As  seen  in  the  figure, 
the  residual  errors  begin  relatively  high,  but  level  off  to  less  than  0.05  after  the  first 
20-30  days.  It  follows  that  the  solutions  obtained  using  iteration  were  not  different 
than  the  linearized  solution.  Table  5.1  shows  that  the  solutions  for  the  linearized 
method,  and  the  iterative  method  with  ^  —  0.05  and  0.02  are  virtually  the  same. 
The  iterative  procedure  requires  more  computations,  as  reflected  by  the  longer  CPU 
times. 

Using  Iteration  with  Larger  Time  Steps  Iteration  can  be  exploited  to  reduce 
simulation  run  times  by  allowing  the  use  of  a  larger  time  step.  Runs  were  made  with 
larger  to  see  if  the  desired  6w:curacy  could  be  obtained  by  using  iterations.  /?  = 
0.02  was  used  for  Af  =  2,  3,  4,  and  5  days.  As  seen  in  Table  5.2,  the  use  of  iteration 
resulted  in  more  acceptable  solutions  for  time  steps  of  2,  3,  and  4  days,  as  compared 
with  the  linearized  solutions  shown  in  Table  4.6.  Any  savings  in  the  larger  time  step 
for  Af  =  3  and  4  days  were  offset  by  the  very  high  number  of  iterations  required.  For 
A/  =  5  days,  even  the  use  of  iteration  did  not  produce  satisfactory  results.  Using 
the  iterative  technique  with  At  —  2  days,  however,  resulted  in  significant  CPU  time 
savings  and  nearly  identical  results  to  the  A<  =  1  day  case.  For  this  ca.se  it  is  clearly 
better  to  use  At  =  2  days  with  iteration,  than  to  use  a  linearized  method  with  At  ~ 
1  day,  since  the  total  run  time  is  approximately  50%  less. 
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Residual  Errortn  Linearized  Solution  for  Baseline  Casa 


Figure  5.1:  Residual  Error  in  Linearized  Solution  of  PDEs,  Baseline  Case 


Table  5.1:  Effect  of  Iteration  on  Baseline  Cane 


TIME 

MEASURE 

Linearii«d 

^  =  1  day 

d  =  0.05 

9  =  0.02 

500  days 

Smmw 

1.74 

1-74 

1.74 

C 

879 

879 

878 

1000  day. 

0.70 

0.70 

0.70 

334 

334 

333 

1500  days 

0.23 

0.23 

0.23 

84 

84 

84 

Totals  for 

Time  to  ilnd  (dayf) 

805 

803 

802 

1500  Days 

Total  /lenitont 

0 

5 

32 

CPU  Tim*  (min) 

72 

98 

100 
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A  Case  Where  an  Iterative  Solution  is  Necessary  Consider  a  case  where  the 
value  of  A'j  was  taken  at  its  minimum  value  of  1  ^  (Table  4.2).  At  this  value,  the 
degradation  occurs  at  a  much  faster  rate  as  seen  by  equations  2.36  -  2.38.  As  seen 
in  Figure  5.2,  the  residual  errors  for  substrate  and  oxygen  are  much  greater  than  for 
the  baseline  case.  The  residual  errors  for  biom^iss  were  extremely  high,  and  were  not 
even  plotted  on  this  graph.  In  this  Ccise,  it  is  obvious  that  a  linearized  solution  is  not 
sufficient. 

Figure  5.3  shows  the  residual  errors  for  the  solutions  obtained  using  iteration  (/?  = 
0.02,  miifr  =  10).  In  this  case  the  error  is  still  very  high  for  the  first  5-10  days,  but 
quickly  levels  off  tr  approximately  5%  for  the  remainder  of  the  1500-day  simulation. 
A  summary  of  the  linearized  versus  the  iterative  solutions  is  given  in  Tabic  5.3. 
Although  this  extreme  example  required  a  relatively  high  number  of  iterations  (and 
long  CPU  time),  the  differences  in  the  predicted  concentrations  were  not  trivial.  In 
fact,  the  very  high  number  of  iterations  was  required  to  improve  the  residual  from 
approximately  0.05  to  0.02.  An  alternative  approach  to  maintaining  accuracy  in  this 
case  of  faster  kinetics  would  be  to  reduce  the  time  step.  As  seen  above,  however,  the 
iterative  method  is  more  efficient  than  the  use  of  a  smaller  time  step. 

Summary  The  use  of  an  iterative  procedure  offers  a  significant  advantage  over  the 
linearized  approach.  In  many  cases  it  permits  the  use  of  a  larger  time  step  than 
would  have  been  otherwise  required.  In  other  cases  where  the  specified  biological 
parameters  result  in  fast  kinetics,  an  iterative  solution  is  required  to  avoid  the  use  of 
a  very  small  time  step.  Thus,  it  is  recommended  as  an  improvement  to  BI02D.  It 
will  be  used  in  all  sensitivity  analysis  runs  completed  for  Chapter  6. 
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Table  5.2:  Using  Iteration  with  Larger  Time  Steps 


Iteration 

=  0.02)  and  At  (days)  of: 

TIME 

MEA,SURE 

1 

2 

3 

4 

5 

500  days 

•S'max 

1.74 

1.73 

1.70 

1.67 

1.47 

Srrxata 

879 

876 

846 

836 

687 

1000  days 

Smax 

0.70 

0.70 

0.67 

0.65 

0.53 

Sma$$ 

334 

331 

314 

304 

233 

1500  days 

^max 

0.23 

0.23 

0.21 

0.20 

0.11 

84 

83 

75 

71 

36 

Totals  for 

Time  io  aind  (days) 

805 

800 

783 

772 

690 

1500  Days 

Total  Iterations 

100 

103 

4500 

2344 

2188 

CPU  Time  (min) 

100 

53 

140 

110 

102 
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Figure  5.2:  Residual  Error  in  Linearized  Solution  of  PDEs,  A' 


5.2  Modified  Kinetics 


The  proposed  changes  to  BI02D  presented  in  this  section  are  the  result  of  discussions 
with  Dr.  James  Gos.  ett,  Cornell  University.  Some  of  the  ideas  that  follow  were  taken 
from  notes  from  Dr.  Gossett’s  class  on  Environmental  Engineering  Processes. 

A  shortcoming  in  BI02D  as  discussed  and  tested  so  far  is  that  the  governing  PDE 
for  oxygen  fails  to  account  for  all  oxygen  sinks.  Specifically,  equation  2.37  neglects 
to  account  for  substrate  consumed  in  synthesis  of  net  biomass  and  O2  consumed  in 
the  aerobic  degradation  of  background  carbon.  Equation  2.37  can  he  rewitten  as: 


Xro  =  BR,F^ 


(5.10) 


where 


= 


f^max 


5 

0 

A', +  5  +  1^ 

[Ko  +  o\ 

(5.11) 


Looking  at  each  additional  O2  sink  individually: 

O2  Consumed  in  Net  Biomass  Synthesis  The  following  equations  represent  the 
production  and  use  of  energy  in  aerobic  respiration  and  synthesis  of  new  biomass, 
assuming  that  the  contaminant  substrate  is  used  as  the  electron  donor  for  both  energy 
production  and  synthesis  of  net  biomass: 

Energy: 


Substrate  +  02—^  CO2  +  H2O  +  energy 

Synthesis: 


(5.12) 


Substrate  +  NH^  +  energy  — >  CzH-j02N  (5.13) 

In  essence,  the  substrate  electron  donor  is  oxidized  and  the  O2  electron  acceptor  is 
reduced.  Energy  is  captured  (e.g.,  as  ATP)  for  use  in  biosynthesis,  which  involves 
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converting  the  substrate  into  cellular  constituants  {€^1^02^).  In  this  case,  when 
the  same  substrate  is  used  for  both  energy  and  synthesis,  a  fraction  of  the  substrate 
is  used  for  synthesis,  while  the  remainder  is  used  for  energy: 


A  +  /«  =  i 


(5.14) 


fa  =  fraction  of  substrate  used  for  synthesis 
fe  =  fraction  of  substrate  used  for  energy 

where  ft  and  /,  depend  upon  energetics  (i.e.,  the  energy  available  from  the  respi¬ 
ration  reaction  compared  to  that  required  by  the  synthesis  reaction).  In  an  activated 
sludge  application,  they  also  depend  upon  the  solids  retention  time  for  the  system. 

In  the  governing  equations  used  in  BI02D,  F,  the  ratio  of  oxygen  to  substrate 
utilized,  assumes  that  all  substrate  is  used  in  an  energy  reaction  (equation  5.12),  and 
that  none  is  used  for  synthesis  (equation  5.13).  They  also  neglect  to  account  for  O2 
used  by  the  microorganisms  as  they  consume  dead  biomass.  An  additional  term  that 
should  be  added  to  (equation  5.10)  is  proposed  to  account  for  these  effects: 


-  1.42 


0 

Ko  +  0 


(5.15) 


where  1.42  is  a  conversion  factor  from  grams  biomass  to  grams  oxygen,  assuming  a 
biomass  composition  of  C5Ht02N.  Thus,  when  there  is  a  net  increase  in  biomass, 
the  (  ]  is  positive,  and  the  amount  of  O2  utilized  will  be  reduced  by  (equation  5.15) 
from  the  amount  otherwise  predicted.  When  there  is  a  net  decreaise  in  biomaiss, 
however,  the  [  ]  is  negative,  and  there  will  be  more  oxygen  consumed  than  in  the 
present  formulation.  The  term  predicts  a  Monod-type  consumption  of  oxygen 

as  the  microorganisms  consume  themselves. 


O2  Consumed  by  Background  Carbon  The  second  oxygen  sink  not  accounted 
for  by  the  governing  PDEs  is  the  O2  consumed  by  the  aerobic  degradation  of  back¬ 
ground  carbon.  Our  simplifying  assumption  about  the  subsurface  microorganisms 
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was  that  they  could  be  represented  by  a  single  facultative  bac.  nai  popi.''.lion  that 
maintained  a  constant  background  concentration  by  aerobically  or  anae-ut-i''.  My  de¬ 
grading  the  available  carbon  sources  in  the  aquifer  (Section  2.3,  Assumpt;-  >’  i).  In  the 
presence  of  a  contaminant  substrate  they  preferentially  consume  it.  P''  '1  ■,  however, 
neglects  to  account  for  any  O2  consumed  by  the  biomass  in  the  aerobic  degradation 
of  this  background  carbon.  In  the  case  where  we  assumed  an  initial,  background 
concentration  of  O2  in  the  aquifer,  as  in  Degradation  Tests  2-4,  ii  follows  that  the 
microorganisms  would  aerobically  degrade  this  carbon.  An  additional  modification  is 
proposed  to  account  for  this; 

+  lA-p^F'CcKc  (5-16) 

ho  +  0 


where  F'  is  the  stoichiometric  ratio  of  oxygen  to  background  carbon  consumed  in  the 
aerobic  degradation  of  background  carbon. 

Thus,  the  modified  PDE  for  oxygen  should  include  the  following: 


^RO  =  I 


HsFB 

Y 


-  1.42 


- 


0 


Ko  +  0 


ThB 


+ 


0 


KoYO 


■'CcA'c} 


(5.17) 


Modification  Testing  The  ultimate  test  of  this  proposed  modification  would  be 
a  comparison  of  predicted  to  observed  concentrations  from  a  field  application.  This, 
however,  is  beyond  the  scope  of  the  present  work.  Tests  were  performed  to  inves¬ 
tigate  how  these  proposed  changes  affect  the  predicted  contaminant  concentrations 
for  Degradation  Test  3  (see  Section  4.4).  This  scenario  involved  the  cleanup  of  a 
aquifer  polluted  by  a  phenol  spill.  The  aquifer  cleanup  was  completed  using  a  single 
injection  well  that  introduced  oxygenated  water  at  a  flow  rate  of  4  The  predicted 
concentrations  for  this  case  were  nearly  the  same  for  BIOPLUME  and  BI02D. 

Figure  5.4  shows  the  predicted  centerline  concentrations  of  phenol  at  500,  1000 
and  1500  days  for  the  original  kinetics  and  the  modified  kinetics.  An  initial  oxygen 
concentration  of  3  ^  was  specified  everywhere  in  the  aquifer,  except  for  anoxic 
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Table  5.4:  Degradation  Test  3:  Modified  vs  Original  Kinetics 


Initial  Aerobic 

Initial  Anoxic 

TIME 

MEASURE 

Original 

ModiHed  1 

1 

Original 

Modified 

500  days 

Smax 

1.74 

1.91 

1.79 

1.91 

^masi 

879 

1402 

1210 

1515 

1000  days 

Smar 

0.70 

1.05 

0.87 

1.06 

Sfnai» 

334 

1212 

803 

1310 

1500  days 

3-nax 

0.23 

0.72 

0.51 

0.73 

Smats 

84 

1161 

574 

1253 

Totals  for 

Time  to  stnd  (days) 

805 

1053 

890 

1070 

conditions  at  the  plume’s  location  (see  Figure  4.5).  Mean  values  for  all  biological 
parameters  were  used  (Table  4.2).  Observe  that  the  differences  in  the  solutions  of  the 
governing  equations  are  magnified  with  time.  This  is  due  in  part  to  the  consumption 
of  the  background  Oj  by  the  facultative  bacteria  when  modified  kinetics  are  used. 
Very  quickly,  the  background  O2  is  consumed  and  anoxic  conditions  prevail.  At  500 
days,  the  initial  O2  is  nearly  gone,  with  the  exception  of  O2  recharged  by  the  injection 
well,  and  O2  recharged  at  the  left  boundary. 

A  more  interesting  comparison  is  seen  in  Figure  5.5.  Again,  the  predicted  cen¬ 
terline  concentrations  of  phenol  at  500,  1000  and  1500  days  are  depicted  for  the  two 
different  kinetic  formulations.  In  this  case,  the  initial  background  O2  condition  was 
taken  as  anoxic  everywhere  i"  the  aquifer.  Here,  the  differences  between  the  predicted 
solutions  were  not  as  great  as  before,  but  still  quite  significant.  The  primary  source 
of  this  difference  is  the  additional  O2  sink  given  by  equation  5.15,  which  accounts  for 
O2  consumed  by  the  consumption  of  dead  biomass. 

Table  5.4  summarizes  the  results  of  these  test  runs. 


Effect  of  Microbial  Kinetics  Model  on  Blo2d  Baseline  Resulls 


Figure  5.4:  Original  vs  Modified  Kinetics  for  Initial  Aerobic  Conditions 


Degradation  Test  3:  Anoxic  Initial  and  Boundary  Cortditlons 


Figure  5.5:  Original  vs.  Modified  Kinetics  for  Initial  Anoxic  Conditions 
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5.3  Summairy 

In  this  chapter  two  modifications  to  the  present  formulation  of  BI02D  were  presented 
and  evaluated: 

1.  An  iterative  procedure  to  solve  the  nonlinear  PDEs 

2.  Modified  kinetics  that  more  realistically  account  for  O2  consumed  in  microbial 
degradation 

The  first  modification  is  highly  recommended  as  a  more  efficient  and  accurate 
means  of  solving  the  problem.  In  some  cases,  where  the  biological  parameters  result 
in  slow  or  moderate  kinetics,  the  iterative  solution  will  permit  more  efficient  solution 
and  shorter  run  times  that  the  linearized  method.  In  Ccises  where  the  kinetics  are 
fast,  the  use  of  an  iterative  solution  is  necessary  to  avoid  the  use  of  an  extremely 
small  time  step. 

The  second  modification  needs  further  testing.  The  inclusion  of  additional  O2 
sinks  due  to  synthesis  of  net  biomass  and  aerobic  degradation  of  background  carbon 
results  in  vastly  different  solutions.  Although  the  modified  kinetics  seem  to  the  author 
as  a  better  representation  of  reality,  model  predictions  must  be  compared  to  observed 
concentrations  at  a  field  site  in  order  judge  this  modification’s  merits.  Until  this  is 
completed,  no  change  is  warranted.  It  is  recommended,  however,  that  the  initial  and 
boundary  conditions  for  O2  in  the  aquifer  be  taken  as  anoxic  everywhere,  unless  there 
is  evidence  to  support  another  assumption.  The  assumption  of  a  constant  aerobic 
initial  condition,  in  the  presence  of  facultative  bacteria,  without  recharge  simply  does 
not  make  sense. 


Chapter  6 

Model  Sensitivity  to  Biological 
Parameters 

6.1  Introduction 

This  chapter  addresses  sensitivity  of  B102D  to  changes  in  the  often  unknown  bio¬ 
logical  parameters  given  in  Table  2.2.  The  purpose  is  to  quantify  the  uncertainty  in 
the  predicted  concentrations  of  pollutant  in  an  aquifer  caused  by  uncertainty  in  the 
estimates  of  these  constants.  Sensitivity  analysis  is  an  essential  step  in  all  applica¬ 
tions  of  groundwater  modeling  [Anderson  and  Woessner,  1992],  but  it  is  especially 
important  in  degradation  modeling  because  of  the  great  number  of  unknown  param¬ 
eters.  As  discussed  in  Chapter  2,  the  predicted  concentrations  are  a  function  of  seven 
physical  parameters  (Table  2.1),  eleven  biological  parameters  (Table  2.2),  and  the 
initial  and  boundary  conditions  for  the  problem.  It  would  be  of  great  help  to  users  of 
degradation  models  to  gain  a  better  understanding  of  which  parameters  are  very  im¬ 
portant  to  know,  eind  which  parameters  are  not  as  important  and  might  be  assumed 
from  literature.  This  will  allow  limited  resources  of  time  and  money  to  be  spent  on 
determining  the  most  important  ones. 


100 


The  basis  for  the  sensitivity  analysis  performed  was  Degradation  Test  3  (see  Sec¬ 
tion  4.4).  This  scenario  involved  the  cleanup  of  a  phenol  polluted  aquifer  using  a 
single  injection  well  that  introduced  oxygenated  water  at  8  Recall  that  for  this 
case  predictions  of  contaminant  concentrations  were  quite  similar  for  BI02D  and 
BIOPLUME.  A  time  step  of  At  =  1  day  was  used.  The  iterative  procedure,  with 
original  kinetics,  presented  in  Chapter  5  was  used  for  all  runs,  with  iteration  param¬ 
eters  of  ^  =  0.02  and  mxitr  =  10. 

Methodology  The  methodology  used  in  this  study  of  model  sensitivity  was  to 
utilize  existing  parameter  ranges  for  phenol  degradation,  observing  predicted  con¬ 
taminant  concentrations  based  on  mean,  meiximum  and  minimum  parameter  values 
from  published  studies.  This  seems  to  be  a  much  more  rational  approach  to  sensi¬ 
tivity  than  merely  varying  each  parameter  by  the  same  percentage,  as  is  often  done. 
By  using  a  realistic  range  the  greater  variablility  of  some  parameters  is  accounted 
for.  For  example,  it  is  much  more  useful  to  know  how  model  predictions  vary  as 
K,  is  changed  from  its  minimum  (1^)  to  its  maximum  value  (266^),  than  it  is  to 
know  how  the  predictions  vary  as  K,  is  perturbed  by  5%  from  its  mean  value.  The 
use  of  the  parameter  extremes  also  allows  for  the  definition  of  bounds  on  the  system 
performance  ("best”  and  "worst”  cases). 

The  ranges  of  the  site-specific  biological  parameters  assumed  for  this  study  are 
given  in  Table  4.2.  These  constants  for  phenol  degradation  were  obtained  from  pub¬ 
lished  laboratory  and  field  studies.  [Ranch  et  al.  ,1983;  Rozich  ct  al.,  1985;  Chang 
and  Rittmann,  1987;  Speitel  et  al.,  1987;  Hobson  and  Mills,  1990;  Lin,  1992;  and 
Wagner,  1992]  Of  the  eleven  parameters,  F,  Cc,  and  fcc  were  considered  as  deter¬ 
ministic.  F  can  be  determined  directly  from  the  balanced  reaction  stoichiometry, 
assuming  complete  mineradization  to  CO2  and  H2O  from  equation  4.4.  Insufficient 
data  were  available  to  quantify  variablity  in  Cc  and  fcc.  and  it  was  assumed  that 


growth  rates  on  background  carbon  were  small  compared  to  growth  on  contaminant 
substrate.  The  remaining  eight  parameters  were  considered  as  unknown  or  stochaistic. 
All  eleven  were  assumed  to  be  constant  in  time  and  space. 

Both  individual  sensitivity  and  combined  sensitivities  were  studied.  The  former 
is  the  more  common  approach  where  model  output  is  observed  as  one  parameter  at  a 
time  is  varied.  In  the  latter  approach  all  parameters  are  varied  simultaneously.  This 
is  more  thorough  and  accounts  for  parameter-parameter  interactions. 

Previous  Work  Taylor  performed  a  sensitivity  analysis  using  probability  density 
functions  for  each  of  the  unknown  parameters.  A  Latin  hypercube  sampling  scheme 
was  used  to  obtain  vectors  of  parameter  distributions.  These  were  input  and  the  vari¬ 
ations  in  model  output  measured.  He  used  Smax  as  the  measure  of  model  sensitivity. 
This  work  differs  from  that  of  Taylor  in  the  following  ways: 

•  Both  individual  and  combined  sensitivities  were  studied 

•  Parameter  extremes  were  used  rather  than  assumed  probability  distributions 

•  Three  measures  of  sensitivity  were  used:  5moi,  Smaas,  and  i,t„j  as  defined  in 
Section  5.1 

•  A  njore  realistic  scenario  was  modeled;  a  larger  aquifer  (80,000  m^)  with  a 
longer  cleanup  time  (1500  days) 

Preliminary  Analysis  A  preliminary,  qualitative  analysis  of  how  the  predicted 
concentrations  should  be  affected  by  changes  in  each  of  these  biological  parameters 
follows: 

1.  F,  the  stoichiometric  ratio  of  oxygen  to  substrate  consumed  in  degradation:  As 
F  increases,  the  amount  of  O2  required  for  degradation  increases  per  unit  of 
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substrate.  Thus,  a  higher  F  should  result  in  iess  degradation  given  a  finite 
amount  of  oxygen.  In  this  work,  however,  F  is  considered  deterministic. 

2.  Umaxi  the  maximum  specific  growth  rate  of  bacteria:  As  pmaz  incre2ises,  the 
biomass  should  grow  more  quickly  and  substrate  consumption  will  procede  at 
a  faster  rate,  resulting  in  more  degradation. 

3.  Kg,  the  Monod  substrate  half-saturation  constant:  As  K,  decreciscs  relative  to  5, 
the  degradation  should  procede  at  a  faster  rate,  resulting  in  more  degradation. 
For  5  »  Kg,  -*  1. 

4.  Ko,  the  Monod  oxygen  half-saturation  constant:  As  Ko  decreases  relative  to  O, 

the  degradation  should  procede  at  a  faster  rate,  resulting  in  more  degradation. 
For  0  »  Ko,  1. 


5.  Ki,  the  substrate  inhibition  constant:  As  A'i  increases  relative  to  5^,  the  bacteria 
are  less  inhibited  by  the  substrate,  and  degradation  should  be  faster.  For  « 

TTi 


6.  Y,  the  yield  coefficient:  As  Y  increases,  more  microbes  are  produced  from  a 
given  amount  of  substrate.  In  addition,  since  Y  is  assumed  constant  for  both 
degradation  of  substrate  and  of  background  carbon,  the  background  biomass 
that  can  be  supported  will  also  increase  with  Y  according  to  equation  4.5. 
Considering  both  effects,  degradation  should  increase  with  higher  values  of  Y. 


7.  rj,  the  endogenous  decay  coefficient:  As  rj  increases,  the  bacteria  decay  at  a 
feister  rate.  In  addition,  an  increase  in  rj,  will  result  in  a  lower  background 
biomass  concentration  according  to  equation  4.5.  The  net  effect  is  unclear; 
however,  it  is  postulated  that  the  amount  of  degradation  will  be  reduced  with 
increased  rj. 
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8.  Rk,  ihe  microbial  retardation  factor:  As  Ri,  increases,  more  biomass  is  adsorbed 
to  the  solid  matrix  of  the  aquifer  and  more  degradation  should  occur.  However, 
according  to  equation  4.5,  increases  in  will  be  offset  by  lower  background 
biomass  concentrations,  and  this  results  in  less  degradation.  The  overall  effect 
is  unknown. 

9.  A'i,  the  Freundlich  isotherm  constant:  As  A'j  increases,  R^  increases,  but  the 
effect  of  increased  Aj  is  unclear. 

10.  iVj,  ihe  Freundlich  isotherm  exponent:  As  N\,  increases,  Aj  increases,  but  the 
effect  of  increased  Aj  is  unclear.  For  N),  ^  1,  A*  =  Ri{B). 

11.  Cc,  the  background  carbon  concentration  in  aquifer:  As  Cc  increases,  the  amount 
of  background  biomass  that  can  be  supported  increases.  Thus  degradation  will 
also  increase.  However,  in  this  study  Cc  is  considered  deterministic. 

12.  kc,  the  first  order  decay  rate  of  background  carbon:  As  kg  increases,  the  amount 
of  background  biomass  that  can  be  supported  increases.  Thu;}  degradation  will 
also  increase,  but  as  with  Ce,  kc  is  considered  deterministic. 

6.2  Sensitivity  to  Uncertainty  in  Individual 
Biological  Parameters 

The  first  part  of  the  sensitivity  study  was  to  determine  the  changes  in  model  predic¬ 
tions  given  the  uncertain  biological  parameters,  considering  one  parameter  at  a  time. 
This  is  the  typical  approach  taken  in  sensitivity  analysis  [Anderson  and  Woessner, 
1992J.  For  the  given  scenario,  simulations  were  performed  with  each  biological  pa¬ 
rameter  at  its  maximum  and  minimum  value,  holding  all  other  f)arameter3  at  their 
mean  values  (Table  4.2).  This  method  has  the  advantage  of  being  straightforward 
and  relatively  easy  to  perform.  For  the  eight  uncertain  parameters,  a  total  of  17  simu¬ 
lation  runs  were  required.  See  Table  D.l  for  a  complete  listing  of  runs  and  parameter 
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combinations.  This  analysis  provided  a  good  initial  estimate  of  what  parameters  were 
most  significant. 

The  results  of  this  analysis  are  summarized  in  Tables  6.1  and  6.2,  and  Figures 
6.1  -  6.3.  The  contaminant  concentrations  predicted  by  BI02D  were  most  sensitive 
to  the  Monod  half  saturation  constant  Kf  The  results  were  moderately  sensitive  to 
fimaxi  '■4,  and  A'o.  As  seen  in  Table  6.2,  the  was  most  sensitive  to  (in  order  of 
increasing  sensitivity)  low  A',,  low  A'o,  low  rj,  high  /iman  high  A',,  and  low  Umax- 
Sensitivity  to  all  other  parameter  extremes  was  less  than  1%  of  for  the  baseline 
case,  and  it  appears  that  at  least  for  this  initial  study,  the  model  results  are  not 
sensitive  to  K,  A',,  A'4,  or  Aj.  The  mean  and  standard  deviation  for  all  measures  was 
computed,  showing  a  re<isonable  match  to  the  baseline  run. 

As  seen  in  equations  2.36  -  2.38,  the  most  sensitive  extremes  (low  A',,  low  A'o, 
low  rj,  high  /imax)  correspond  to  faster  microbial  kinetics.  In  these  runs,  the  greatest 
number  of  iterations  in  the  solution  of  the  nonlinear  PDFs  were  required  for  these 
cases  with  fast  kinetics.  Thus  it  appears  that  the  degree  of  nonlinearity  may  be  pro¬ 
portional  to  the  rate  of  degradation.  It  particular,  notice  in  Table  6.1  the  extremely 
high  number  of  iterations  required  for  run  #5,  corresponding  to  A’,  =  1  In  many 
of  the  1500  times  steps  the  iterative  procedure  did  not  converge  to  a  2%  tolerance  in 
the  10  iterations  allowed.  See  Section  5.1  for  a  more  complete  analysis  of  this  case. 

Comparing  these  results  to  what  was  expected  (preliminary  an2dysis)  it  is  seen 
that  changes  in  K„  fimax,  and  rj  caused  model  predictions  as  anticipated.  Changes 
in  concentrations  resulting  from  Y,  Ki,  A'j,  eind  jVj  were  not  evident,  leading  to 
the  initial  conclusion  that  the  model  predictions  are  insensitive  to  these  parameters. 
Model  results  from  Kg  extremes  were  not  as  predicted.  For  A'o  =0.1  less 
degradation  was  predicted  than  for  the  baseline  case,  however,  the  model  predictions 
for  Ao  =  2  ^  were  quite  close  to  the  baseline  run.  This  was  not  well  understood 
until  further  analysis  of  sensitivity  changes  with  time  was  completed  (see  below). 
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Indtvlduat  Sensitivity  Results 


Figure  6.1:  Individual  Sensitivity:  versus  Run  (Runs  Defined  in  Table  D.l) 

It  is  intersesting  to  note  that  of  the  seventeen  runs,  three  predicted  more  degra¬ 
dation  than  BIOPLUME  [K,  =  1  =  15.36  3^,  and  rj  =  0.001  3^^).  One 

must  conclude  that  either  BI02D  is  overpredicting  degradation  or  that  BIOPLUME 
is  not  accurately  depicting  instant^eous  kinetics,  due  to  numerical  errors. 

The  sensitivity  of  the  model  output  to  A',,  timaxy  and  Kg  is  further  illustrated  in 
Figures  6.4  -  6.7.  These  plots  depict  centerline  concentration  profiles  at  500  and  1500 
days  for  BI02D  as  each  parameter  was  varied  from  its  minimum  to  its  maximum 
vadue.  Each  plot  also  shows  how  these  extremes  compare  with  BIOPLUME.  It  is 
evident  from  these  plots  that  using  Smax  as  the  only  index  of  sensitivity  m.ay  not 
be  sufficient.  It  is  important  to  measure  sensitivity  by  both  5fna*  and  Smaix-  For 
example,  for  run  5^4  where  K,  —  266  ^  (Figure  6.4),  5m«i  at  500  days  is  only 
2%  greater  than  the  baseline  case;  but  Smass  is  56%  greater.  The  estimate  of  total 
contaminant  mass  was  more  sensitive  to  A',  than  was  the  estimate  of  Smax- 
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Table  6.1:  Sensitivity  To  Changes  in  Individual  Biological  Parameters.  (Runs  Defined 
in  Table  D.l) 


PARAMETER 

500  days 

1000  days 

1500  days 

Iterations 

RUN 

VALUE 

Smax 

5 

SmatB 

^Mtnd 

required 

1 

Baieline 

1.74 

879 

0.70 

334 

0.23 

84 

805 

100 

2 

=  15.38 

1.65 

754 

0.65 

262 

0.19 

49 

765 

59 

3 

Umax  —  1.97 

1.77 

1155 

0.72 

499 

0.26 

183 

816 

30 

4 

K,  =  266 

1.78 

1369 

0.74 

661 

0.30 

292 

822 

29 

5 

>: 

II 

1.34 

477 

0.43 

104 

0.02 

1 

636 

13118 

6 

Ki  =  1463 

1.74 

877 

0.70 

333 

0.23 

84 

802 

32 

7 

Ki  =  23 

1.74 

879 

0.70 

333 

0.23 

84 

803 

32 

8 

Y  =  1.92 

1.74 

877 

0.70 

333 

0.23 

84 

802 

32 

9 

y  =  0.50 

1.74 

877 

0.70 

333 

0.23 

84 

802 

32 

10 

r»  =  0.10 

1,76 

999 

0.71 

403 

0.24 

125 

810 

44 

11 

r»  =  0.001 

1.62 

687 

0.64 

228 

0.20 

35 

759 

25 

12 

1.74 

928 

0.70 

362 

0.23 

102 

800 

31 

13 

K,  =  0.1 

1.78 

964 

0.83 

459 

0.39 

172 

867 

50 

14 

=  22.97 

1.75 

888 

0.71 

338 

0.23 

86 

808 

32 

{Rt  =150) 

15 

Kt  =  7.55 

1.74 

874 

0.70 

331 

0.22 

83 

800 

32 

{Rt  =50) 

18 

Ni  =  1.1 

1.75 

883 

0.70 

336 

0.23 

85 

805 

0 

(R*  =  99B^  +  1) 

17 

Nk  =  0.9 

1.75 

883 

0.70 

338 

0.23 

85 

805 

0 

{Rt  =  99B*  +  1) 

MEAN 

1.71 

897 

0.69 

352 

0.23 

101 

794 

STND  DEV 

0.10 

185 

0.08 

117 

0.07 

65 

47 

BIOPLUME 

1.72 

756 

0.70 

243 

0.22 

41 

795 

Smu  at  SOO  Day*  (mg/L| 
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Table  6.2:  Most  Sensitive  Parameters  in  Individual  Study 


RUN 

PARAMETER 

VALUE 

^»tnd 

(days) 

%  from 

BASELINE 

5 

A'.  =  1 

636 

-21.0 

13 

=  0.1 

867 

+7.7 

11 

r»  =  0.001 

759 

-5.7 

2 

^‘max  =  15.36 

765 

-5.0 

4 

K,  =  266 

822 

-(■2.1 

3 

/^max  “  1*97 

816 

-1-1.4 

Individual  Sansitivtty  Results 


i 


Figure  6.2:  Individual  Sensitivity:  Smax  versus  Run  (Runs  Defined  in  Table  D.l) 
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Individual  Sensitivity  Results 


Figure  6.3:  Individual  Sensitivity;  Smaas  versus  Run  (Runs  Defined  in  Table  D.l) 


Substrate  Half  Saturation  Constant  Extremes 


Figure  6.4:  Sensitivity  Analysis:  Effect  of  K,  Extremes  (Runs  4,5) 
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Sensitivity  Over  Time  P  ensitivity  index  as  described  by  Fjeld  et  al.  [1989] 
was  used  in  order  to  further  study  the  most  sensitive  parameters.  This  method 
involves  using  normalized  measures  of  Smax  and  Smass  that  emphasize  changes  on 
an  absolute  scale.  Large  percentage  variations  are  not  significant  when  S  <C  S,fndi 
whereas  small  percentage  changes  may  be  quite  significant  when  S  S,t„d.  The  goal 
of  this  analysis  was  to  investigate  how  sensitivity  changes  with  time.  Let  T]t  and  71 
be  defined  aa  follows; 

cextreme _ cmean 

Vt  = - Q -  (6.1) 

“Bind 

ceTtreme _ cme&n 

7,  =  ZBiSh - (6.2) 

SmMB,ave 

where  5ma**’"*  is  the  maximum  contaminant  concentration  in  the  aquifer  at  time  <, 
found  by  a  simulation  with  the  parameter  at  its  mciximum  or  minimum  value. 
is  the  maximum  contaminamt  concentration  in  the  aquifer  at  time  t,  found  by  a  sim¬ 
ulation  with  all  parameters  at  their  mean  values  (baseline  case).  Similarly, 
is  the  total  contaminant  mass  in  the  aquifer  at  time  t,  found  by  a  simulation  with  the 
p2urameter  at  its  maximum  or  minimum  value.  is  the  total  contaminant  mass 

in  the  eiquifer  at  time  t,  found  by  a  simulation  with  all  parameters  at  their  mean 
values  (baseline  case).  It  is  now  possible  to  plot  r/j  and  7t  versus  time  to  see  when 
the  model  is  especially  sensitive.  In  this  study  Smax  Wcis  normalized  with  respect  to 
S$indi  and  Smass  was  normalized  with  respect  to  the  average  mass  in  the  aquifer,  but 
in  other  applications  different  selections  may  be  more  appropriate. 

Figures  6.8  -  6.15  show  how  the  sensitivity  indices  r/j  and  71  vary  with  time  for 
K„  A'o,  and  rj  extremes.  A  negative  index  indicates  that  the  model  prediction 

using  the  parameter  extreme  is  smaller  than  the  prediction  made  with  all  parameters 
at  their  mean  values. 
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Summary  Observations 

•  For  rjt,  the  normalized  sensitivity  index  of  Smaz  in  equation  6.1,  the  model 
predictions  were  not  sensitive  to  low  pmai,  high  K,,  high  rj,  nor  high  Ko- 
This  WM  not  true,  however,  for  7,  the  normalized  sensitivity  index  of  Smass 
in  equation  6.2.  7i  was  sensitive  to  low  and  high  values  of  fimax,  high  r>. 
This  demonstrated  that  there  are  some  significant  differences  in  the  sensitivities 
depending  on  the  performance  measure  of  interest. 

•  r}t  sensitivity  to  high  Umax,  low  A',,  low  rj,  and  low  Ko  was  very  high  initially 
and  decreased  with  time.  For  high  fimax  and  low  rj,  r/j  -+  0  with  increasing 
time.  This  is  not  true  for  low  A'j,  whose  sensitivity  is  still  significant  at  1500 
days. 

•  Sensitivity  of  rjt  and  7<  to  low  Ko  demonstrated  a  reversal  at  approximately  400 
days.  For  the  first  400  days,  <  0  and  7j  <  0,  which  is  consistant  with  what 
was  expected.  After  400  days,  the  effect  of  low  Ko  was  that  rjt  >  0  and  7j  >  0. 
This  may  be  due  to  the  fact  that  the  value  of  A'o  is  really  only  significant 
relative  to  O.  Thus,  in  the  early  stages  of  the  cleanup  A'o  <  0,  but  as  the 
background  oxygen  is  utilized,  Ko  >  O.  This  results  in  a  kinetic  limitation 
and  less  degradation.  This  crossover  effect  is  further  illustrated  in  Figure  6.16. 
Observe  that  at  200  days  the  centerline  profiles  of  phenol  were  as  expected 
with  the  concentration  predicted  with  A'o  =  2  greater  than  the  concentration 
predicted  with  Ao  =  0.1.  At  400  days  the  concentrations  predicted  for  the  two 
Ko  extremes  were  approximately  equal.  At  600  days,  however,  the  effect  of  Kg 
has  reversed  and  the  concentrations  predicted  with  Kg  =  0.1  were  greater  than 
with  A'o  =  2. 


S«isifcitir^lnd«.Ela  _  d  ^  Conorntraljon  (tngrt.) 
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Oxygon  Half  Saturation  Constant  Extremes 


!  6.7:  Sensitivity  Analysis:  Effect  of  A'o  Extremes  (Runs  12,13) 


Sanaitlvity  of  Smax  vs  Time 


Time  (days) 


Figure  6.8:  Sensitivity  Index  rft  (Equation  6.1)  Versus  Time  for  fimaz  (Runs  2,3) 
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Sensitivity  of  Total  Mass  vs  Tima 


Figure  6.9:  Sensitivity  Index  7t  (Equation  6.2)  Versus  Time  for  Umax  (Runs  2,3) 


Sensittvity  of  Smax  vs  Tima 


Figure  6.10:  Sensitivity  Index  T]t  (Equation  6.1)  Versus  Time  for  K,  (Runs  4,5) 
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Sanstttviry  of  Total  Mass  vs  Time 


Figure  6.11:  Sensitivity  Index  7j  (Equation  6.2)  Versus  Time  for  K,  (Runs  4,5) 


Sensitivity  of  Sm«x  vs  Time 


Figure  6.12;  Sensitivity  Index  r/j  (Equation  6.1)  Versus  Time  for  rj  (Runs  10,11) 
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Sensttlvtty  of  Smass  vs  Tima 


Figure  6.13:  Sensitivity  Index  ft  (Equation  6.2)  Versus  Time  for  rj  (Runs  10,11) 


SansItMty  of  Smaj<  vs  Tima 


Figure  6.14;  Sensitivity  Index  vt  (Equation  6.1)  Versus  Time  for  Ko  (Runs  12,13) 


117 


SansHivfty  of  Smass  vs  Time 


Figure  6.15:  Sensitivity  Index  74  (Equation  6.2)  Versus  Time  for  Ko  (Runs  12,13) 


Oxygen  Half-Saturation  Constant  Extremes 


Figure  6.16:  Combined  Sensitivity:  Changing  Effects  of  Kg 


118 


6.3  Combined  Sensitivity  to  Uncertainty  in 
Biological  Parameters 

The  analysis  of  BI02D’s  sensitivity  to  uncertainty  in  individual  hhlogica.]  parameters 
yielded  some  valuable  information,  but  it  w^is  incomplete.  The  approach  of  varying 
one  parameter  at  a  time  assumed  no  parameter- parameter  interactions.  But  this  is 
not  realistic  since  there  is  uncertainty  in  all  parameters.  A  more  complete  method  is 
to  consider  joint  or  combined  sensitivity. 

2*  Factorial  Design  The  method  selected  for  the  study  of  joint  effects  was  bor¬ 
rowed  from  a  field  of  Operations  Research  called  Experimental  Design  and  Opti¬ 
mization,  as  presented  by  Law  and  Kelton  [1991],  They  describe  a  computationally 
economical  strategy  to  study  the  respones  of  a  simulation  model  given  possible  de¬ 
cision  levels  or  factors.  The  method  is  called  2*  factorial  design.  In  the  context 
presented  by  Law  and  Kelton,  this  is  a  decision  making  tool  used  to  achieve  optimal 
performance  in  a  given  system.  Here  managers  or  analysts  are  using  a  simulation 
model  to  assess  the  effects  of  two  or  more  possible  changes  or  improvements,  and 
want  to  consider  decision  interactions. 

2*  factorial  design  provides  a  systematic  way  of  measuring  model  sensitivity  where 
these  interactions  are  important.  The  method  is  conceptually  straightforward;  for 
each  decision  level  or  factor  k,  two  possible  values  or  design  points  are  selected,  and 
a  total  of  2*  simulation  runs  are  performed.  Law  and  Kelton  also  define  quantitative 
measures  of  the  factor  main  effects  and  of  two  and  three-way  interactive  effects. 

The  method  of  2*  factorial  design  was  directly  applied  to  BI02D  and  groundwater 
biodegradation  modeling.  The  uncertain  biological  parameters  were  taken  as  the 
design  factors  (while  their  /ilues  are  not  controlled  or  designed,  selection  of  reasonable 
values  allow  for  application  of  the  method).  For  each  parameter  the  two  possible 
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values  were  taken  as  the  maximum  and  minimums  for  phenol  (Table  4.2).  The  use 
of  2*  factorial  design  proved  to  be  a  very  useful  way  of  measuring  and  comparing  the 
combined  sensitivity  of  BI02D  to  multiple  changes  in  the  biological  parameters. 

Assumptions  As  in  the  study  of  individual  sensitivity,  the  parameters  F,  Cc,  and 
kc  were  considered  deterministic.  Based  on  the  results  of  the  individual  sensitivity 
analysis,  it  was  further  assumed  that  microbial  adsorption  could  be  modeled  with 
a  linear  isotherm  (jVj  =  1),  and  that  inhibition  effects  were  insignificant  (A'i  ^ 
S)  so  that  A',  was  taken  cis  constant  at  its  mean  value.  These  assumptions  were 
made  in  order  to  reduce  the  number  of  required  simulations  from  2®  to  2®,  and 
were  supported  by  results  from  the  earlier  study.  Thus,  a  total  of  64  runs  were 
simulated,  corresponding  to  the  64  different  combinations  of  parameters.  The  runs 
were  organized  in  a  systematic  manner  so  that  run  pairs  were  easily  identified.  A 
run  pair  corresponds  to  two  runs  where  one  parameter  is  varied  from  its  minimum 
to  its  maximum  value,  while  the  other  five  parameters  were  held  constant.  Thus, 
there  were  a  total  of  32  run  pairs  for  each  parameter.  A  complete  listing  of  runs  and 
parameter  combinations  is  found  in  Tables  D.2  -  D.3. 

Results  The  results  of  the  combined  sensitivity  runs  are  presented  in  Tables  6.3 
-  6.5,  and  Figures  6.17  -  6.25.  Tables  6.2  -  6.4  show  the  model  predictions  of  Smax 
Juid  Smasa  for  500,  1000  and  1500  days  and  the  total  time  required  to  achieve  a  1  ^ 
standard  {tgind)  everywhere  in  the  aquifer. 

Model  predictions  of  t,i„d  and  of  Smax  and  Smasa  for  500  days  are  presented  as 
scatter  plots  in  Figures  6.17  -  6.19.  A  close  look  at  these  plots  reveals  some  definite 
patterns  that  correspond  to  specific  parameter  combinations.  For  example  observe 
Figure  6.17.  Here,  the  first  32  runs  were  made  with  fimax  =  1-97  and  the  last  32 
with  fi-nax  =  15.36  In  runs  33-64,  there  is  a  definite  downward  shift  of  the  even 
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Combined  Sensitivity  Results 
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Figure  6.17;  Combined  Sensitivity:  ittni  versus  Run  (Runs  Defined  in  Tables  D.2  - 
D.3) 

numbered  runs  (corresponding  to  K,  —  266  seen  above  the  baseline,  whereas 
the  odd  numbered  runs  (A',  =  1  do  not  shift  as  dramatically  with  the  change  in 
A'max*  Thus,  low  Kg  appears  to  dominate  over  Umax  because  the  results  for  the  the 
cdd  numbered  runs  do  not  change  much  with  changes  in  ftmaz-  For  high  Kg,  however, 
the  effect  of  changes  in  pmai  were  observed  as  more  dramatic.  This  supports  earlier 
observations  that  the  model  is  most  sensitive  to  Kg. 

The  low  sensitivity  of  results  to  A'j  can  be  seen  by  compauing  runs  1-16  with  runs 
17-32,  and  runs  33-48  with  49-64.  Observe  in  Figures  6.17  -  6.19  that  there  is  almost 
no  noticeable  difference  between  the  results  from  runs  pairs  1,17;  2,18;  3,19;  etc.  The 
pattern  in  each  quadrant  repeats  itself  with  only  three  exceptions  (pairs  43,59;  45,61; 
and  47,63).  Thus,  it  appears  that  the  model  is  insensitive  to  Aj,  especially  at  low 
values  of  Umax- 
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Table  6.3:  Combined  Sensitivity  To  Changes  in  Biological  Parameters  for  Runs  1-30 
(Runs  Defined  in  Table  D.2  -  D.3) 


I 


Table  6.4:  Combined  Sensitivity  To  Changes  in  Biological  Parameters  for  Runs  31-60 
(Runs  Defined  in  Table  D.2  -  D.3) 


500  days 

tas  Smi 


Jmaat 


632.27 

2377.50 

125.72 

781.30 

222.92 

786.11 
217.32 

761.30 

215.11 

786.11 

122.30 
1193.09 

0.00 

1431.98 

452.05 

1192.47 

504.40 

1431.97 

153.66 

754.24 

205.02 

776.92 
125.14 
754.24 
208.49 

778.92 
58.62 

1190.81 

856.21 

1431.67 


lOCX)  days 

Sttiat  SmasM 

0.544 

180.61 

1.139 

2026.49 

0.000 

0.00 

0.737 

361.46 

0.026 

0.99 

0.713 

279.27 

0.003 

0.09 

0.737 

361  .<46 

0.019 

0.65 

0.713 

279.27 

0.000 

0.00 

0.822 

540.06 

0.000 

Q.QQ 

0.751 

n4.i6 

0.190 

34.32 

0.822 

539.64 

0.454 

114.52 

0.751 

724.15 

0.000 

0.00 

0.998 

500.00 

0.008 

0.20 

0.705 

273.94 

0.000 

0.00 

0.998 

500.00 

0.015 

0.42 

0.705 

273.94 

0.000 

0.00 

0.821 

538.63 

0.852 

294.90 

0.751 

724.01 

1500  days 
Smax  Smai 


Iterations 

retfaired 


13492 

48 

1000 

26 

1437 

20 

908 

26 

1092 

20 

8424 

111 

3847 

48 

11884 

114 

13496 

48 

1052 

28 

1278 

21 

1070 

26 

1272 

21 

7125 
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Table  6.5:  Combined  Sensitivity  To  Changes  in  Biological  Parameters  for  Runs  61-64 
(Runs  Defined  in  Table  D.2  -  D.3) 


RUN 

500  days 

1000  days 

^max  Srna$t 

1500  days 

Smax  Sma»i 

^itnd 

Iteration* 

required 

61 

0.632 

132.37 

0.000 

0.00 

0.000 

0.00 

407 

8489 

S3 

1.851 

1189.92 

0.820 

538.01 

0.341 

198.41 

874 

91 

63 

0.939 

288.11 

0.197 

28.08 

0.000 

0.00 

473 

11599 

64 

1.784 

1431.62 

0.7S1 

723.98 

0.491 

500.00 

829 

48 

2.5 


Comb<n«d  Sensitivity  Results 
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Figure  6.18;  Combined  Sensitivity;  Smaz  versus  Run  (Runs  Defined  in  Tables  D.2  - 
D.3) 


Figure  6.19:  Combined  Sensitivity:  Smass  versus  Run  (Runs  Defined  in  Tables  D.2 
D.3) 
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Main  Effects  This  qualitative  analysis,  however,  is  difficult,  and  can  be  somewhat 
misleading.  Law  and  Kelton  present  a  way  of  quantifying  these  results  as  follows.  Let 
efc  be  the  main  effect  of  parameter  k.  It  is  defined  as  the  average  change  in  model 
prediction  due  to  changing  parameter  k  from  its  minimum  to  its  maximum  value. 
This  average  is  taken  over  all  combinations  of  the  other  parameters  in  the  problem, 
ejfc  can  alternately  be  defined  as  the  difference  between  the  model  predictions  with 
k  at  its  maximum  value  and  the  model  predictions  with  k  at  its  minimum  value. 
It  is  possible  to  define  and  ej,""'*  using  Smax,  Smaas  and  ijtnd  as  the 

respective  performance  measures.  For  and  the  main  effect  is  averaged 

over  the  entire  1500  day  simulation.  For  example,  the  main  effect  of  A',  on  is 


defined  as: 


J.ind  _  (^Stnd  -  ~  ^ainj)  +  •  •  •  + 


'A".  - 


where  1,2, 3,  •  ■  •  ,64  refer  to  the  simulation  run  number  (Appendix  D).  Thus,  for  A', 
the  run  pairs  were  (2, 1),  (4,3),  •  •  • ,  (64,63).  Similarly,  the  main  effect  of  ftmax  on 


is  defined  as: 


El  500 
<=i 


The  main  effect  of  Y  on  Smasa  is  defined  as: 


E1500  ^rr 

<=1  - - 


Using  these  and  similar  definitions,  main  effects  of  the  six  uncertain  parameters  on 
the  three  performance  measures  were  calculated.  As  presented  in  Table  6.6,  a  positive 
effect  can  be  explained  as  an  increase  in  performance  measure  as  the  parameter  is 
changed  from  its  minimum  to  its  meiximum  value.  For  example,  as  Kg  is  changed 
from  1  ^  to  266  the  main  effect  was  an  increase  in  predicted  cleanup  time  of 
383.7  days  and  aji  increase  in  the  average  Smax  in  the  aquifer  of  0.6023 
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Table  6.6:  Combined  Sensitivity:  Main  Effects  (Given  By  Example  in  Equations  6.3 
-  6.5) 


PARAMETER 

k 

Main  Effect: 

^Smax 

ef-*“  [g) 

ej*'"''  [days) 

Kg 

0.6023 

848.2 

384 

f^max 

-0.0682 

-334.8 

-115 

rb 

0.2992 

459.9 

81 

Kg 

0.2158 

103.7 

74 

Y 

-0.0259 

13.3 

13 

Kb 

0.0600 

1.1 

7 

The  maun  effects  in  Table  6.6  indicate  that  the  predicted  concentrations  of  con¬ 
taminant  are  most  sensitive  to  A%,  the  Monod  half  saturation  constamt.  The  study 
results  aJso  indicate  that  the  model  is  relatively  insensitive  to  changes  in  Y  and  A'j. 
For  fimaxi  ^01  J«id  rj  the  results  are  moderately  sensitive,  and  the  order  of  sensitivity 
changes  with  the  measure  of  interest.  For  t,tndj  ^he  order  of  sensitivity  was  A',,  Umax, 
rj,  and  Ko\  for  Smax  the  order  was  AT,,  r^,  Kg,  and  Umax,  and  for  5Tna»*  the  order 
was  A",,  rj,  Umax,  and  Kg. 


Bar  Graphs  An  easier  way  of  depicting  these  results  is  to  plot  the  difference  be¬ 
tween  the  predicted  performance  measures  with  the  parameter  at  its  minimum  value 
and  with  the  p.:.rameter  at  its  maocimum  value,  while  all  others  were  held  constant  at 
one  particular  combination.  This  was  done  for  each  parameter  pair  using  as  the 
performance  measure.  The  results  are  plotted  as  bar  charts  in  Figures  6.20  -  6.25. 
In  each,  the  straight  line  indicates  ejj.""'*,  the  main  effect,  which  is  the  average  of  all 
bars  shown  (and  is  given  by  example  in  Equation  6.3). 
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Several  observations  are  made  about  these  plots. 

•  The  first  32  runs  (16  pairs),  corresponding  to  Umax  =  1-97  produced  pat¬ 
terns  that  were  nearly  symmetric.  See  bars  to  the  left  of  the  vertical  line  in 
Figures  6.20,  and  6.22  -  6.25.  In  these  runs,  there  was  almost  no  sensitivity  to 
changes  in  Y  or  A'j. 

•  In  the  l<ist  32  runs,  however,  where  /imaz  =  15.36  3^  much  more  variation  in 
the  results  was  seen.  See  bars  to  the  right  of  the  vertical  line  in  Figures  6.20, 

.  and  6.22  -  6.25.  This  was  particularly  noticeable  for  odd  numbered  runs,  which 
correspond  to  A',  =  1  This  was  probably  due  to  the  combined  effects  of 
a  high  Umax  and  a  low  AT,;  both  extremes  cause  faster  degradation  (equations 
2.36  -  2.38). 

•  Runs  43,  45,  47  and  59  produced  unexpected  results.  This  was  particularly 
noticeable  as  the  extremely  high  or  low  bars  in  Figures  6.20  -  6.25.  These  pair 
numbers  are  highlighted  on  the  bcir  charts.  The  common  trait  in  each  of  these 
cases  was  a  high  fiman  a  low  Ks,  and  a  high  rj.  Again,  this  corresponds  to 
faster  microbial  kinetics. 

•  All  plots  were  made  with  the  same  scale  on  the  Y-axis,  so  that  the  relative 
sensitivity  can  be  seen  qualitatively  by  comparison  of  the  bar  graphs. 

•  Ko  effects  were  drastically  different  for  run  pairs  4,2;  8,5;  20,18,;  and  24,22 
where  ^imax  =  l-S’f  3^)  Al,  =  266  ^  and  rj  =  .001  (Figure  6.22).  In  these 
run  pairs,  the  effect  was  negative,  whereas  in  all  other  runs  (of  the  first  32)  the 
effect  was  positive.  This  effect  wais  also  seen  in  pairs  44,42;  48,46;  60,58;  and 
64,62  where  fimax  =  15.36  3^,  K,  =  266  ^  and  rj  =  .001  3^.  Thus,  for  low 
rj,  the  main  effect  of  Kg  is  opposite  in  direction  to  A',  {  —  Kg  at  high  A”,,  and 
-l-Ao  at  low  A',).  But  at  high  rj,  the  main  effect  of  Kg  is  positive  regardless  of 

A,. 
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Comb<n»d  Ssnsitivity:  Ks 


Figure  6.20:  Combined  Sensitivity:  Differences  in  t,tnd  for  A',  Pairs,  As  Defined  in 
Table  D.4 

Interactive  Effects  The  main  effects  presented  above  calculate  the  average  change 
of  model  output  due  to  changes  in  uncertain  biological  parameters.  This  average  is 
taken  over  all  the  various  combinations  of  1:  —  1  other  parameters.  It  is  possible  to 
measure  the  degree  of  interaction  between  parameter  ki  and  k2  by  the  use  of  two-way 
interaction  effect^  ^*1*3  defined  as  one  half  the  difference  between  the  average 

effect  of  parameter  ki  when  parauneter  ^2  is  at  its  maximum  level  and  the  average 
effect  of  parameter  ki  when  parameter  k2  is  at  its  minimum  level.  Another  way  of 
defining  the  two-way  interaction  effect  is  the  difference  between  the  average  effects  of 
parameters  /ti  and  k2  when  both  are  at  the  same  extreme  (maximum  or  minimum), 
and  the  average  effects  when  the  are  at  opposite  extremes.  The  two-way  effects  are 
symmetric  (cfcjtj  =  For  example,  the  two-way  interaction  effect  of  K,  and 
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Combmod  Sansitivrty:  Ko 


Figure  6.22:  Combined  Sensitivity:  Differences  in  tstnd  for  Ko  Pairs,  As  Defined 
Table  D.4 
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Umax  using  as  the  measure  of  performance  is: 


g**,<n<J  _  2. 

K$iimax  O 


^33 


64 


«63 


16 


(6.6) 


i^ltnd  ~  ^\ind)  + - ^  i^stnd  ~ 


16 


where  1, 2, 3,  ■  ■  • ,  64  refer  to  the  simulation  run  number  (Appendix  D).  The  interaction 
effect  of  Ks  and  A'o  using  Smax  the  measure  of  performance  is: 


-K.Ko,  - 


i^maxt  -  ■^mor,)  -f - h  {Smax,  ~  -^mazt) 

16 

i^max,  ~  ^Laxt)  ^ - i^maxt  ~  •^mai,) 


16 


(6.7) 
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Smax 

hiha 


EISOO  *5mar 
t=l  ^K.Kof 


1500 


(6.8) 


Several  of  the  various  two-way  interactive  effects  are  shown  in  Table  6.7.  For  the 
performance  measures  of  Smax  aJid  Ks  —  Kn  was  the  strongest  two-way  effect. 
For  the  performance  measures  of  Smasa,  however,  the  effects  of  Ks-f^mcx  and  A',- 
rj  were  the  most  important.  The  magnitude  of  the  two-way  effects  can  be  directly 
compared  to  the  main  effects  calculated  previously  (Table  6.6).  Observe  that  the 
two-way  effects  are  in  general  smaller  than  the  main  effects. 

Three-way  effects  were  also  computed  in  a  maniier  similar  to  the  two-way  effects. 
The  three-way  effects  are  also  symmetrical  (etjtjtj  =  ^k^kikj  —  ^k^kiks^  ®tc.  )  Three- 
way  effects  for  A",  —  Umax  —  ^uid  A,  —  A©  —  rj  were  computed  and  are  shown  in 
Table  6.7.  They  are  in  general  smaller  than  the  main  effects  and  two-way  effects. 
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Table  6.7:  Combined  Sensitivity:  Interactive  Effects  (Given  by  Example  in  Equations 
6.6  -  6.8) 


PARAMETERS 

Interactive  Effects: 

*1  -  *2 

(5) 

“  Pmar 

-0.138 

-264.5 

-33 

“  Pmajr 

-0.019 

-15.6 

-42 

K.-Ko 

-0.188 

-41.0 

-72 

K.-Y 

0.026 

-13.4 

-13 

K.-n 

-0.135 

336.9 

7 

K.-K, 

-0.063 

-4.5 

-10 

-0.081 

-196.4 

-34 

h\  -  rj 

0.079 

100.1 

43 

-0.038 

-85.6 

2 

“  rj 

-0.007 

29.9 

22 

6.4  Conclusions 

The  sensitivity  analysis  performed  in  this  chapter  yielded  many  new  insights  into 
BI02D,  and  how  its  predictions  of  contauninant  concentrations  vary  with  uncertainty 
in  biological  parameters.  These  observations  are  summarized  below: 

•  Individual  sensitivity  alone  was  insufScient  in  trying  to  assess  the  effects  of 
uncertainty  in  BI02D.  Parameter-parameter  interactions  were  found  to  be  sig¬ 
nificant.  As  seen  in  Table  6.8,  the  effect  of  considering  individual  sensitivity 
alone  would  have  been  to  misjudge  Kg  as  the  second  most  sensitive  parameter 
and  misjudge  the  direction  of  the  sensitivity  of  A'o. 

•  A  high  coefficent  of  variation  for  an  uncertain  parameter  does  not  necessarily 
lead  to  high  sensitivity.  For  example,  the  coefficient  of  variation  for  A',  was 
1.14,  second  only  to  K,,  but  BI02D  was  not  sensitive  to  this  parameter. 
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•  The  contaminant  concentrations  predicted  by  BI02D  were  most  sensitive  to 
chajiges  in  A’j,  the  Monod  half  saturation  constant.  Model  sensitivity  to  Kg 
is  undoubtedly  due  in  part  to  the  variablity  of  this  but  this  was  true  for  both 
combined  and  individual  sensitivity,  and  for  all  three  performance  measures. 
This  observation  is  also  consistant  with  what  Tay/or  found  in  his  study.  Thus,  in 
application  of  BI02D,  it  is  recommended  that  resources  in  parameter  estimation 
be  spent  on  determining  K,  ever  any  otner  uncertain  pcirameters. 

•  The  contaminant  concentrations  predicted  by  BI02D  were  not  sensitive  to 
changes  in  A'j,  Y ,  nor  Ki.  Thus  it  is  recommended  that  resources  in 
parameter  estimation  not  be  spent  determining  these  parameters.  Values  taken 
from  literature  (for  the  specific  contaminant  of  interest)  should  be  sufficient  for 
these. 

•  The  sensitivity  of  the  predicted  contaminant  concentrations  to  Ko  was  quite 
variable,  depending  upon  time,  and  the  values  of  the  other  constants. 

•  Results  of  the  sensitivity  analysis  differed  with  the  performance  measure  of 
interest.  For  example,  the  importance  of  ptmax  is  greater  when  one  is  interested 
in  estimating  the  time  to  meet  a  1  ppm  cleanup  standard,  than  knowing  Ko 
or  rj.  However,  if  the  tot2J  contaminant  mass  is  chief  concern,  then  it  is  more 
important  to  know  rj  than  Umax  or  Kg. 

•  Combined  effects  of  K^-Ko  were  most  significant  of  the  parameter  pairs  when 
^tind  or  Smax  was  the  performance  measure.  However,  in  estimating  the  total 
contauninant  mass,  the  combined  effects  of  K,-  Umax  and  A'j-rj  were  the  most 
important. 
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Table  6.8:  Comparison  of  Individual  and  Combined  Sensitivity  for  (Given  by 
Example  in  Equation  6,3) 


PARAMETER 

CoefBcient  of 

Variation 

Main  Effect 

on  t,tnd  (days): 

Individual 

Sensitivity 

Combined 

Sensitivity 

1.32 

186 

384 

f^max 

.55 

-51 

-115 

n 

1 

51 

81 

Ko 

1 

-67 

74 

Y 

.25 

0 

13 

Ki, 

.50 

6 

7 

Nb 

.10 

0 

- 

K„ 

Nb 

Ki 


.10 

1.14 


-1 


Chapter  7 


Summary 

7.1  Conclusions 

In  this  work  an  in-situ  bioremediation  model  was  tested  and  studied  in  great  detail. 
Specifically,  BI02D  was  validated  in  its  ability  to  model  flow  and  transport  using 
the  IGWMC’s  Level  1  Testing  Protocol.  It  was  then  compared  to  the  better  known 
and  proven  model  BIOPLUME  in  its  ability  to  model  biodegradation.  Two  potential 
improvements  to  BI02D  were  presented  and  evaluated.  Fineilly,  BI02D  was  studied 
to  determine  which  input  biological  parameters  were  most  important  in  determining 
predicted  concentrations.  Based  on  this  work,  the  following  conclusions  are  made: 

•  BI02D  did  an  acceptable  job  modeling  flow  and  transport  for  the  five  IGWMC 
test  cases  applied  to  it  in  Chapter  3.  While  there  are  certainly  more  sophis¬ 
ticated  models  available,  BI02D  is  a  good,  relatively  straightforward  model 
that  allows  the  additional  complications  of  degradation  to  be  included  without 
making  it  too  cumbersome. 

•  As  with  many  numerical  codes,  BI02D  is  prone  to  numerical  errors  if  the  ac¬ 
curacy  criteria  given  in  equations  2.33  -  2.35  are  not  satisfied. 
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•  BI02D  did  a  good  job  of  modeling  biodegradation  in  the  four  test  cases  eval¬ 
uated.  It  predicted  slightly  less  degradation  than  BIOPLUME  in  three  of  the 
four  degradation  tests  performed  in  Chapter  4,  as  Wcis  expected. 

•  For  the  degradation  test  involving  an  injection-pumping  well  doublet,  BI02D 
actually  outperformed  BIOPLUME  based  on  numerical  dispersion  and  mass 
balance  errors. 

•  The  use  of  an  iterative  procedure  to  solve  the  nonlinear  PDEs  is  recommended 
over  a  linearized  solution.  This  method  was  proven  as  more  efficient  in  Chapter 
5. 

•  Individual  sensitivity  analysis  alone  was  proven  as  insufficient  for  BI02D.  An 
more  complete  and  accurate  method  was  presented  that  included  the  effects  of 
parameter-parameter  interactions. 

•  BI02D  was  most  sensitive  to  the  Monod  half  saturation  constant,  A'j.  It  was 
not  sensitive  to  A'j,  Nf,,  Y,  nor  A',-.  Sensitivity  to  rj,  A'o,  and  pmoi  was  moderate 
and  depended  to  some  extent  upon  the  performance  measure  of  interest.  There¬ 
fore,  in  application  of  BI02D,  it  is  recommended  that  resources  in  parameter 
estimation  be  spent  on  determining  A",  over  any  other  uncertain  parameters. 

7.2  Recommendations  for  Further  Research 

The  following  areas  need  further  attention  and  are  recommended  for  future  research: 

•  The  application  of  BI02D  to  a  real  site  is  the  only  way  to  truly  validate  this 
model.  It  will  also  permit  further  evaluation  of  the  choice  of  kinetics  discussed 
in  Chapter  5.  The  validation  of  BI02D  at  one  of  the  sites  that  BIOPLUME 
was  applied  is  a  reasonable  pos.sibility.  although  these  two  cases  involve  natural 
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degradation  (no  pumping)  only.  This  is  the  best  way  to  gain  confidence  in 
BI02D’s  ability  to  model  the  complexities  of  biodegradation. 

•  Further  sensitivity  analysis  should  be  performed  using  the  combined  approach 
presented  in  Chapter  6.  It  is  recommended  that  an  alternate  selection  of  the 
two  values  used  in  2*  factorial  design  be  tested.  In  addition  to  the  minimum 
and  maximum,  another  measure  related  to  the  parameters’  real  distributions 
such  as  ±one  standard  deviation  should  be  tested.  This  will  serve  to  verify  the 
results  of  Chapter  6  and  to  avoid  the  problems  associated  with  combinations 
of  extremes. 

•  Sensitivity  analysis  of  the  optimization  model  would  be  an  interesting  and  im¬ 
portant  work.  It  would  be  of  great  value  to  learn  how  the  optimal  cleanup 
strategy  chosen  changes  with  uncertainty  in  the  parameters. 


Appendix  A 


Definition  of  Variables  Used 


VARIABLE 

DESCRIPTION 

UNITS 

T 

transmissivity  tensor 

[^] 

K 

permeability  tensor 

[^1 

h 

vertically  averaged  hydraulic  head 

[m] 

S 

storativity 

[1] 

b 

saturated  aquifer  thickness 

[m] 

n 

porosity 

[m] 

C 

concentration  of  conservative  tracer 

5 

concentration  of  substrate  (degradable  contaminent) 

[71 

0 

concentration  of  oxygen 

[^1 

B 

concentration  of  microbial  biomass 

m 

B, 

background  concentration  of  biomass 

m 

D 

dispersion  tensor 

[^] 

Do 

molecular  dispersion 

[3^] 
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OL 

longitudinal  dispersivity 

[m] 

QX 

transverse  dispersivity 

[m] 

V 

Darcy  velocity 

^RS 

rate  of  substrate  biodegradation 

ll^l 

^RO 

rate  of  oxygen  consumption 

l-A] 

^RB 

rate  of  biomass  growth  (+)  or  decay  (-) 

Qw 

injection  (+)  or  extraction  (-)  rate  in  BI02D 

[3^] 

injection  (-)  or  extraction  (+)  rate  in  BIOPLUME 

[^1 

Su, 

substrate  concentration  in  water  source/sink 

[7] 

Ow 

oxygen  concentration  in  water  source/sink 

[^1 

By, 

biomass  concentration  in  water  source/sink 

[^1 

R, 

substrate  retardation  factor  = 

[1] 

Ro 

oxygen  retardation  factor  =  1 

[11 

Ri 

biomass  retardation  factor  =  ^1  + 

[1] 

Pb 

soil  bulk  density  (taken  as  =  2.65) 

Kd 

substrate  partitioning  coefficient 

[^1 

Kb 

Freundlich  isotherm  partitioning  coefficient 

for  bacteria  adsoption 

[^1 

Ni 

Freundlich  isotherm  exponent  for  bacteria  adsoption 

[1] 

F  ratio  of  oxygen  to  substrate  used  in  degradation 

F'  ratio  of  oxygen  to  background  carbon  used  in  degradation 

Y  yield  coefficient 

P-mux  maximum  specific  growth  rate  for  bacteria 
fia  bacterial  growth  rate 

rj  bacterial  decay  rate 

Ks  Monod  half-saturation  coefficient  for  substrate 

Kg  Monod  half-saturation  coefficient  for  oxygen 

Ki  Inhibition  coefficient  for  substrate 

Cc  natural  organic  carbon  in  aquifer 

kc  first  order  decay  rate  of  organic  carbon 
Sij  Dirac  delta  function  evaluated  at  (ij) 

A'o  vertical  conductivity 

Wi,  Wj  finite  element  basis  functions 


[^] 

[^] 

(^1 

(il 

[afc] 

[2^] 

m 

m 

[^] 

1^1 

[sfc] 

(11 


n  unit  outward  normal  vector  to 


finite  element  boundary 


[1] 


^t)  ny  finite  element  directional  cosines  between  unit 
normal  to  boundary  and  x,y  coordinate  axes 
nfe  number  of  finite  elements  in  aquifer  mesh 
nds  number  of  nodes  in  aquifer  mesh 
nbw  band  width  of  finite  element  matrices 


[1] 

[1] 

[11 

[1] 
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Pe 

Peclet  Number 

[1] 

Cr 

Courant  Number 

[1] 

^max 

maximum  substrate  concentration  in  aquifer 

m 

^maaa 

total  substrate  mass  in  aquifer 

b] 

^itnd 

time  to  meet  a  1[^]  standard 

[days] 

0 

Convergence  criterion  for  iterative  procedure 

[1] 

6a 

residual  error  in  substrate  equation 

[1] 

So 

residual  error  in  oxygen  equation 

[1] 

Si, 

residual  error  in  biom^iss  equation 

[1] 

mxitr 

maximum  number  of  iterations  in  iterative  procedure 

[1] 

V 

normalized  sensitivity  measure  of  Smaz 

[1] 

7 

normalized  sensitivity  measure  of  Smaaa 

[1] 

e* 

main  effect  of  parameter  k 

variable 

Appendix  B 

BI02D  Finite  Element  Equations 


This  appendix  details  the  finite  element  equations  used  in  the  simulation  model 
BI02D.  The  equations  are  based  on  an  unpublished  manuscript  by  Taylor  [1991], 
and  Taylor  [1993].  All  mat: ices  and  vectors  are  defined  in  Appendix  A. 

The  hydraulic  head  equation  derived  from  equation  3.1  is: 


[Am  =  +  {F,}  (B.l) 

where  Qu,,  is  the  pumping  rate  vector  as  a  function  of  time. 

The  elemental  coefficient  matrices  are  defined  as  follows: 

{nK  =  (/r  “i  •  t  k,  ■  d/j  (B.3) 

ini  =  14 . eri  (B.4) 


where  u’,  is  the  weighting  function  at  node  i,  and  riy  are  directional  cosines  between 
the  unit  normal  to  the  boundary  and  the  x,  y  coordinate  axes,  e\  are  binary  vectors 
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denoting  the  location  of  the  i**  well  in  the  head  vector,  F  is  the  surface  of  the  element, 
and  hj  is  the  specified  boundary  head.  All  other  terms  are  defined  as  in  Appendix 
A. 

The  substrate  (contaminant)  concentration  equation  derived  from  equation  3.2 
using  finite  elements  and  a  variably  weighted  backward  difference  scheme  is 

([M.l  +  M6  [[yv]  +  +  Er=i  QwjPc]'])  =  (8.5) 

([M,]  -  A<(i  -  9)  [[yv]  +  {sy  + 

LT=i  QwJPcYS^  +  Ate{F,y^^  +  Ai(i  - 

where  {5}*  and  are  substrate  concentrations  vectors  at  the  beginning  and 

end  of  the  current  time  step  respectively.  At  is  the  time  step  ,  and  6  is  the  time 
weighting  or  implicitness  factor  (0  <  ^  <  1).  Equation  B.5  is  nonlinear  because  of 
[/?s]{5}.  The  matrices  are  defined  as  above  and 


[yv]«.  =  + 

dx  dy] 


(B.6) 


[M,\-j  =  \^j  J  {bnR,)wiWjdx  dy 


(8.7) 


= 


u 


bnwiWkBlRf, 

k=l 


t  n  Mmai 


1 


A',  +  5i  +  ^ 


Oi 


Wjdxdy 


(8.8) 


=  (^(6A>- V5)-nu-,d/)' 


(8.9) 
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[PcY^eu-el  (B.IO) 

where  fz  is  a  normal  vector  to  the  boundary,  Ci^c  are  binary  vectors  denoting  the 
location  of  the  well  in  the  concentration  vector. 

Similar  to  equation  B.5,  the  oxygen  concentration  equation  derived  from  equation 
3.3  is: 


{[Mo]  +  Ate  [[iV]  +  [/?'+»]  +  E.’l,  {0}'+'  =  (B.ll) 

{[Mo]  -At{l-  9)  [[iV]  +  [fi‘]  +  EIli  {0}'  + 

HZx  +  A<(1  -  e){Foy 

where  {0}*  and  {O}*"*"’  are  oxygen  concentrations  vectors  at  the  beginning  and  end 
of  the  current  time  step  respectively.  Equation  B.ll  is  nonlinear  because  of  [i?o]{0}. 
Additional  matrices  not  defined  above  are: 

["•It  =  [//  (i>n)  uZjWjdx  dyj  (B.12) 


’  r  n  1  ® 

/  j  YlhnwiWkB[Ri,^^^  - ^ ^  ■■■  wjdxdy  (B.13) 


{For,  =  (^J^{bD-VO)-nw,diy 


(B.14) 


Finally,  the  biomass  equation  derived  from  equation  3.4  is: 


+  (1  -  e][M,YAt9  [(.V]  +  +  E;^,  Qu-,,[Fc]’])  =  (B.15) 

{9[M,Y^^  +  (1  -  0)[A/i]‘  +  Ai{l  -  9)  [[iVj  +  [Rl]  +  Er=i  Q^.,[Pc]‘])  {BY  + 
E;1,  Qu,JPc]'Bo,  +  A<0{FJ‘+’  +  A/(l  -  9){FiY 
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where  {B}*  and  {B}*"*"*  are  biomass  concentrations  vectors  at  the  beginning  and  end 
of  the  current  time  step  respectively.  Equation  B.15  is  nonlinear  because  of  [Mft]{5} 
whenever  a  nonlinear  isotherm  is  used.  Additional  matrices  not  defined  above  are: 


[A/j] =  J  J  (bnRi,)  w,Wjdx  dy 


(B.16) 


II 


bnWtWkRh 

k=\ 


f^max 

Y 


Si 


A',  +  5j  +  f 


oi 


Ko  +  Oi 


Wjdxdy 


(B.17) 


{CiK"  Vfl)  •  nwidl  )'-[//  (hnCY kcWi)  dx  dyj 


(B.18) 


Appendix  C 


Validation  Testing:  Problem 
Statements  and  Analytical 
Solutions 

C.l  Validation  Test  1:  Transport  in  a 
Semi-Infinite  Column 

Problem  Statement  One-dimensional  advective-dispersive  transport  of  a  conser¬ 
vative  contaminant  through  a  semi-infinite  porous  medium  is  given  by  the  following 
equation  [Bear,  1979]: 


„  dC 

^  dx^ 


(C.l) 


where  D  is  the  coefficient  of  longitudinad  dispersion,  R,  is  the  retardation  factor,  and 
V  is  the  velocity  in  the  x-direction.  The  initial  and  boundary  conditions  are: 


C(i,0)  =  0 


(C.2) 
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150 

(C.3) 
(C.4) 

where  Co  is  the  constant  concentration  at  the  inlet  boundary. 


C(0,<)  =  Co 
C(oo,  <)  =  0 


Analytical  Solution  Ogata  and  Banks  [1961]  give  the  anlytical  solution  of  the 
problem  as 


C_ 

Co 


1 

2  exp 


Vi 

2D 


—  Vx  Rgnx  —  Vt  Vx  ,  R,nx+Vt 

exp  er/c - +  exp— rrer/c - 

^  2D  ^  2y/RTDi  ^  2D  ^  2^/R]Wt 


(C.5) 


Input  Specifications  The  physical  parameters  used  in  Validation  Test  1  are  shown 
in  Table  C.l  and  Figure  4.1.  The  grid  used  in  this  problem  consisted  of  40  elements 
with  nodal  spacing  of  10  meters  in  the  x-direction.  Thus,  the  overall  length  of  the 
column  was  400  meters.  The  simulation  was  conducted  for  50  days  at  an  initial 
time  step  of  2.5  days.  Other  time  steps  were  considered  in  order  to  maintain  the 
implicitness  factor  at  1.0.  As  shown  in  Table  C.2,  the  simulation  took  1.30  seconds 
on  an  IBM  RS/6000  (Model  570)  workstation. 
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Table  C.l:  Values  of  Physical  Parameters  Used  in  Validation  Test  1 


Parameter 

Value 

Darcy  Velocity,  V 

Porosity,  n 

0.25 

Longitudinal  Dispersivity,  qi 

5  m 

Concentration  at  Boundary,  Co 

1.0  ^ 

Ax 

10.0  m 

Rs  for  Case  1 

1.0 

Rs  for  Case  2 

2.0 

At  (initially) 

2.5  days 

Implicitness  Factor,  & 

variable 

Peclet  number  (Peg) 

2 

Courant  number  (Cr) 

0.4 

Table  C.2:  Validation  Test  1  Timing  Results 


CPU  Time  (sec) 

No.  of  Elements  No.  of  Nodes 

Per  time  step  For  20  time  steps 

40  82  1  0.065  1.3 
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C.2  Validation  Test  2:  Transport  From  A 
Continuous  Point  Source 

Problem  Statement  This  test  involves  the  two-dimensional  dispersion  of  a  con¬ 
servative  solute  injected  from  a  point  source  in  uniform  areal  groundwater  flow.  It  is 
assumed  that  the  injection  rate  is  so  small  that  the  natural  velocity  of  groundwater 
is  unchanged.  For  transport  of  a  conservative  species,  the  governing  equation  is 

„  d'^c  „  ,ac  ^dc  ^ 

Dtt  ^^2  +  Dyy  Qc  (C.6) 

where  Dxz  and  Dyy  are  the  coefficients  of  hydrodynamic  dispersion  in  the  x  and  y 
direction.  Qc  is  the  mass  injection  rate  of  solute  per  unit  volume  of  aquifer  [wl- 
is  also  assumed  that  mechanical  dispersion  dominates  over  molecular  diffusion.  Thus, 
the  coefficients  become 


Dxz  =  a  lV 

Dyy  =aTV 


where  ai  and  aj  are  the  coefficients  of  longitudinal  and  transverse  dispersivity. 
The  initial  and  boundary  conditions  for  this  problem  are: 

C(x,j/,0)  =  0  (C 

<5c(x,t/,0  =  (C 


C(±oo,  ±oo,  f)  =  0 


(C.IO) 


where  Q  is  the  volumetric  injectun  rate  of  fluid  per  unit  aquifer  thickness  -y  ,  Cq 
is  the  concentration  of  the  injected  fluid  J-pj  and  S  is  the  Dirac  delta  function 
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Analytical  Solution  Following  the  procedure  described  in  Bear  [1979],  the  general 
solution  may  be  written  as 


where 


QCpexp^ 

AieJo  XX  ^yy 


(“•5) 


(C.ll) 


(C.12) 


r  =  .  h;2  ^ 


R,nr^ 

«  =  — - 

ADxxt 

^exp  f-  f^  +  7^ 


(C.13) 


(C.14) 


(C.15) 


W{u,  is  commonly  referred  to  as  a  Hantush  [1956]  well  function  for  the  problem 
of  transient  flow  to  a  well  in  an  infinite  leaky  aquifer.  Thus,  ast— ♦oo,oru— ►Oa 
steady-state  condition  will  be  reached.  This  results  in  a  balance  between  the  rate  of 
solute  injection  and  dispersion.  Equation  C.ll  becomes 


c  =  (^) 

AieJDx.Dyy  \B] 


(C.16) 


where  A'o(^)  is  the  modified  Besssel  Function  of  the  second  kind,  zero  order. 


Input  Specific  ations  Values  of  the  physical  parameters  used  in  Validation  Test  2 
are  shown  in  Table  C.3  and  Figure  4.5.  As  discussed  in  Huyakorn  et  al.  [1984]  the 
data  represents  a  simulation  of  the  hexavalent  chromium  contaminant  reported  by 
Perlmutter  and  Lieber  [1970]  and  Wilson  and  Miller  [1978].  This  problem  considered 
three  different  d-icre"' izations  of  the  aquifer:  fine,  medium  and  coarse. 

The  problem  is  axisymetrical  with  respect  to  the  x-axis,  so  it  was  possible  to 
model  the  upper  b  l;'of  the  aquifer  only,  saving  on  memory  requirements  and  speeding 


..'•■■■■  ,.  ■  ■  '  t'  .’  ••,;.•'•••  ^  'V:'-  ■  ■•' 
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Table  C.5:  Values  of  Physical  Parameters  Used  in  Validation  Test  2 


Parameter 

Value 

Aquifer  Thickness,  b 

33.5  m 

Darcy  Velocity,  V 

.161  f 

Porosity,  n 

0.35 

Longitudinal  Dispersivity,  a/, 

21.3  m 

Transverse  Dispersivity,  aj' 

4.3  m 

Contaminent  mass  flux,  QCo 

R. 

1.0 

At  (initially) 

100  days 

Implicitness  Factor,  9 

variable 

Selected  Grid 

Ax 

Ay 

P.r 

Pc, 

Cr 

Fine 

10  m 

10  m 

.47 

2.32 

1.61 

Medium 

60  m 

15  m 

2.8 

3.5 

.27 

Course 

60  m 

30  m 

2.8 

7.0 

.27 

calculations.  The  nodal  spacing,  and  number  of  elements  and  nodes  for  each  case 
are  shown  in  Table  C.4.  The  simulation  took  86.8,  21.2,  and  5.8  seconds  on  an  IBM 
RS/6000  workstation  for  the  fine,  medium  and  coarse  grids  respectively. 
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Table  C.4:  Validation  Test  2  Timing  Results 


Selected  Grid 

No.  of  Elements 

No.  of  Nodes 

CPU  Time  (sec) 

Per  time  step  For  20  time  steps 

Fine 

3000 

3171 

4.34 

86.8 

Medium 

800 

861 

1.06 

21.2 

Coarse 

400 

451 

0.29 

5.8 

C.3  Validation  Test  3:  Transport  Of  A  Solute 
Slug  In  A  Uniform  Flow  Field 

Problem  Statement  This  test  involves  the  two-dimensional  dispersion  of  a  con¬ 
servative  solute  slug  injected  in  uniform,  isotropic  groundwater  flow.  It  is  assumed 
that  the  injection  doesn’t  change  the  natural  groundwater  velocity.  For  transport  of 
a  conservative  species,  the  governing  equation  is 


„  d'^C  „  d'^C 

Dxx  k  T  d" 


Bc  _  p  ac  _ 


(C.17) 


5x2  ■  ^yvdy2 

where  Dxx  and  Dyy  are  the  coefficients  of  hydrodynamic  dispersion  in  the  x  and  y 
direction.  Qc  is  the  mass  injection  rate  of  solute  per  unit  volume  of  aquifer  [;^]-  If 
is  also  assumed  that  mechanical  dispersion  dominates  over  molecular  diffusion,  and 
the  coefficients  are 


Dxx  =  olV  (C.18) 

Dyy  =  OtV 


where  ai  and  or  are  the  coefficients  of  longitudinaland  transverse  dispersivity. 
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The  initial  and  boundary  conditions  for  this  problem  are: 

C(i,y,0)  =  0  (C.19) 

=  (C.20) 

n 

C(ioo,  ±oo,  t )  =  0  (C.21) 

(C.22) 

where  m  is  the  solute  ma^s  injected  per  unit  aquifer  thickness  and  6  is  the  Dirac 
delta  function  [L~'^]. 

Analytical  Solution  The  general  analytical  solution  was  presented  by  Sauty  [1980], 
and  using  present  notation: 


i'p 


c«(a,4)  =  -^f^exp 


+  t 


n 


4t', 


RMAX 


R  J 


where 


,  -r.2 


x%  + 


aiQT 

tji  -  - 


+ 


Vr 


ain 

^'rmax  =  -  2 


XR=  { 


X 


(C.23) 

(C.24) 

(C.25) 

(C.26) 

(C.27) 


I  >  0, 

|i|  +  ^  I  <  0. 

The  concentration  C  is  calculated  as  a  product  of  the  dimensionless  concentration 
Cr  and  the  peak  concentration  C^fAX'- 

M 


CmaxI^r^vr)  = 


where 


f{xR,yR)  = 


1 


XRiRMAX 


exp  — 


4-Knaiy/aiaT 
^fi(l  -  tRMAx)^ 


f{^R,yR) 


‘iiRMAX 


+ 


Vr 


‘^^R^RMAX 


(C.28) 


(C.29) 


157 


Table  C.5:  Values  of  Physical  Parameters  Used  in  Validation  Test  3 


Parameter 

Value 

Aquifer  Thickness,  h 

1.0  m 

Darcy  Velocity,  V 

r>  T7l 
^  1 

Porosity,  n 

0.35 

Longitudin.il  Dispersivity,  ai 

4  m 

Transverse  Dispersivity,  aj 

1  m 

Contaminent  mass  per  unit  thickness  of  aquifer 

3500  ^ 

Rs 

1.0 

Az 

5  m 

Ay 

5  m 

At  (initially) 

1  day 

Implicitness  Factor,  0 

variable 

Peclet  number  (Pe,) 

1.25 

Peclet  number  {Pey) 

5.00 

Courant  number  (Cr) 

0.40 

with 


iRMAX  = 


\xr/  \xrJ 


2_ 

XR 


(C.30) 


Input  Specifications  Values  of  the  physical  parameters  used  in  Validation  Test 
3  are  shown  in  Table  C.5  and  Figure  4.8.  A  rectangular  finite  element  grid  was 
used  in  the  simulation  with  Ai  =  Ay  =  5m.  The  problem  is  axisymetrical  with 
respect  to  the  x-axis,  so  it  was  possible  to  model  the  upper  half  of  the  aquifer  only, 
reducing  memory  requirements  and  speeding  calculations.  As  shown  in  Table  C.6, 
the  simulation  took  19.6  seconds  on  an  IBM  RS/6000  workstation. 
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No.  of  Elements  No.  of  Nodes 

CPU  Time  (sec) 

Per  time  step  For  15  time  steps 

960  1029 

1.31  19.65 

C.4  Validation  Test  4:  Transport  In  A  Radial 
Flow  Field 


Problem  Statement  This  validation  test  involves  the  two-dimensional  dispersion 
of  a  solute  injected  from  a  fully  penetrating  well  in  a  confined  aquifer  (Figure  4.12). 
It  is  assumed  that  the  rate  of  injection  is  constant  and  that  the  background  aquifer 
i^ow  is  negligible.  The  resulting  flow  field  is  assumed  radial  and  steady-state.  The 
following  equation  desribes  the  problem: 


dc 

^  dr^  ^  dr  dt 


(C.31) 


where  r  is  the  radial  distance  from  the  well  and  V  is  the  Darcy  velocity  given  by 


(C.32) 


where  Q  is  the  injection  rate  of  the  fluid  and  h  is  the  aquifer  thickness  [L]. 
The  initial  and  boundary  conditions  for  this  problem  are: 


(C.33) 

(C.34) 

(C.35) 

(C.36) 


C(r,0)  =  0 

C(oo,  t)  =  0 
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Analytical  Solution  The  general  analytical  solution  was  derived  by  Ogata  [1958] 
and  presented  by  Huyakom  et  al.  [1984].  Using  present  notation: 


C 

Co 


where 


'r  - 

too  exp 

.  2ai  . 

Jo  V 

M(u)  dv  (C.37) 


a  = 


a  — 


o')  -h  Y}[a>) 

G  = 

Q 

lithn 

^  3 

2  1 

-  & 

Z-^JaiG 

1/2 

^  3 

2  ( 

Zy/aiG 


(C.38) 

(C.39) 

(C.40) 

(C.41) 


and  Ji  and  are  Bessel  functions  of  order  3,  of  the  first  and  second  kinds  respec¬ 
tively. 

Equation  A. 35  is  too  complicated  to  evaluate  analytically.  Huyakorn  et  al.  [1984] 
recommend  using  Hoopes  and  Harleman’s  [1966]  approximate  solution: 


a = 


't-Gi 


air^i 


(C.42) 


where  G  is  defined  cis  above. 


Input  Specifications  Values  of  the  physical  parameters  used  in  Validation  Test 
4  are  shown  in  Table  C.7.  The  problem  is  axisymetrical,  so  that  it  was  possible  to 
model  one  quarter  of  the  aquifer  only,  saving  considerably  on  memory  requirements 
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Table  C.7:  Values  of  Physical  Parameters  Used  in  Validation  Test  4 


Parameter 

Value 

Well  Discharge,  Q 

25  ^ 

day 

Initial  Concentration,  Co 

Aquifer  Thickness,  b 

10  m 

Darcy  Velocity,  V 

Of 

Porosity,  n 

0.25 

Longitudinal  Dispersivity, 

0.3  m 

Transverse  Dispersivity,  aj 

0.0  m 

R, 

1.0 

Ax 

1  m 

Ay 

1  m 

At  (initially) 

1  day 

Implicitness  Factor,  9 

variable 

Peclet  number  (Pe^) 

3.33 

Courant  number  (Cr) 

0.40 

and  speeding  calculations.  A  rectangular  finite  element  grid  was  used  in  the  sim¬ 
ulation  with  Ax  =  Ay  =  Im.  A  total  of  400  finite  elements  and  441  nodes  were 
used.  As  shown  in  Table  C.8,  the  simulation  took  19.6  seconds  on  an  IBM  RS/6000 
workstation. 
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Table  C.8:  Validation  Test  4  Timing  Results 


No.  of  Elements  No.  of  Nodes 

CPU  Time  (sec) 

Per  time  step  For  40  time  steps 

400  441 

0.49  19.60 

C.5  Validation  Test  5:  Transport  in  a 

Non-Uniform  Flow  Field  Caused  by  an 
Injection-Extraction  Well  Doublet 


Problem  Statement  The  final  validation  test  used  to  evaluate  BI02D  concerns 
solute  transport  between  a  pair  of  recharging  and  discharging  wells  operating  at  a 
constant  flow  rate.  Both  wells  fully  penetrate  a  uniform  thickness  confined  aquifer 
that  is  assumed  as  infinite,  homogeneous  and  isotropic  (Figu.''e  3.14).  The  flow  field 
is  assumed  as  steady-state.  For  transport  of  a  conservative  species,  the  governing 
equation  is 


dx  **  dx  dy  dy  dx  dy  dy  dx  *  ^  dy 


(C.43) 


where  Dix,Dxx,  and  Dxx  are  the  three  components  of  the  hydrodynamic  dispersion 
tensor,  and  I4  and  Vy  are  the  components  of  the  Darcy  velocity  in  the  x  and  y 
directions,  respectively.  For  the  well  doublet,  the  values  of  14  and  I4  are  given  by 


I  -  lO 

I  -f  xo 

(C.44) 

'  2x6 

_(x  -  xo)2  -1- 

(x  -f  xo)^  + 

1 

1  1 

(C.45) 

2x6 

(x  -  xo)2  -f-  y2 

(x  +  xo)2  -f  y^j 

where  Q  is  the  well  flow  rate  [^],  b  is  the  aquifer  thickness  and  xq  is  half  the  well, 
spacing. 

The  initial  and  boundary  conditions  associated  with  equation  C.43  are: 
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C(x,y,Q)  =  0  (C.46) 

C{-xo,y,0)  =  Co  (C.47) 

Analytical  Solutions  For  the  most  general  case  involving  the  combined  influences 
of  advection  and  dispersion,  an  analytical  solution  does  not  exist  [Huyakorn  et  ai, 
1984J-  For  a  more  limited  case  of  pure  advection,  a  solution  was  developed  by  Hoopes 
and  Harleman  [1967]  and  Gringarten  and  Sauty  [1975].  For  a  second  limiting  case  of 
advection  and  longitudinal  dispersion  only,  an  approximate  analytical  solution  was 
developed  by  Hoopes  and  Harleman  [1967]  and  Grove  and  Beetem  [1971],  Both  are 
presented  in  [Huyakorn  et  ai,  1984]- 

For  the  first  case  (pure  advection),  a  semi-analytical  solution  was  developed  and 
programmed  by  Javandel  et  al.  [1984]-  The  model,  called  RESSQ,  uses  the  complex 
velocity  potential  to  estimate  the  concentration  distribution  in  the  aquifer.  Validation 
Test  5  utilized  RESSQ,  which  assumes  a  confined  aquifer  of  uniform  thickness,  infinite 
in  extent,  that  is  homogeneous  and  isotropic. 

The  technique  used  by  Javandel  et  al.  was  to  first  identify  simple  flow  components, 
such  as  sinks  and  sources.  Then  the  overall  complex  velocity  potential  of  the  system 
is  obtained  by  combining  the  expressions  for  each  individual  component.  The  velocity 
field  is  then  identified  by  taking  the  derivative  of  the  velocity  potential.  Locations  of 
contaminant  fronts  and  flow  patterns  are  estimated  at  various  values  of  time.  Finally, 
stream  function  of  the  system  is  used  to  calculate  the  time  variation  of  the  rate  at 
which  a  contaminant  reaches  any  desired  point  [El-Kadi,  1988]. 
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Input  Specifications  Values  of  the  physical  parameters  used  in  Validation  Test 
5  are  given  in  Table  C.8.  The  problem  is  axisymetrical  with  respect  to  the  x-axis, 
so  that  it  was  possible  to  model  one  half  of  the  aquifer  only,  saving  considerably  on 
memory  requirements  and  speeding  calculations.  A  rectangular  finite  element  grid 
was  used  in  the  simulation  with  Ax  =  Ay  =  5  m. 


Table  C.9;  Physical  Parameters  For  Validation  Test  5 


Parameter 

Value 

Injection  and  Extraction  Rates,  Q 

±2.0  ^ 

Injected  Concentration,  Cq 

100  ^ 

Darcy  Velocity,  V 

0.015  = 

Porosity,  n 

0.29 

R, 

1.0 

Aquifer  Thickness,  b 

1.0  m 

Ax 

5.0  m 

At/ 

5.0  m 

At 

1.0  day 

Implicitness  Factor,  9 

1.0 

Case  1 

Dispersivities,  ai  =  ar 

.01  m 

Peclet  number  (Pe,  =  Pcy) 

500 

Case  2 

Dispersivities,  orr,  =  ot 

.1  m 

Peclet  number  (Pe,  =  Pe,) 

50 

Case  3 

Dispersivities,  ai  =  or 

1  m 

Peclet  number  (Pcr  =  Pcy) 

5 

Case  4 

Dispersivities,  ai  =  or 

5  m 

Peclet  number  (Per  =  Pcy) 

1 

Case  5 

Longitudinal  Dispersivity,  Oi 

9.1  m 

Tranverse  Dispersivity,  ox 

1.8  m 

Peclet  number  (Pcr) 

.55 

Peclet  number  (Pe,) 

2.78 

Appendix  D 

Sensitivity  Analysis  Run 
Summaries 

This  appendix  details  the  sensitivity  analysis  runs  discussed  in  Chapter  6.  The 
biological  parameters  are  defined  in  Table  2.2,  and  the  specific  parameter  ranges 
considered  for  phenol  are  found  in  Table  4.2. 

Table  D.l  contains  the  specific  values  of  the  biological  parameters  used  for  the 
individual  sensitivity  analysis  presented  in  Section  6.2.  A  total  of  17  runs  were  made 
in  the  individual  sensitivity  analysis.  The  results  of  these  runs  are  found  in  Tables 
6.1  -  6.2,  and  Figures  6.2  -  6.16. 

Tables  D.2  -  D.3  contain  the  specific  values  of  the  biological  parameters  used  for 
the  combined  sensitivity  analysis  presented  in  Section  6.3  A  total  of  64  runs  were 
made  in  the  individual  sensitivity  analysis.  The  results  of  these  runs  are  found  in 
Tables  6.3  -  6.5,  and  Figures  6.17  -  6.19. 
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Table  D.4  contains  the  specific  run  numbers  that  comprise  parameter  run  pairs 
used  in  Section  6.3  and  Figures  6.20  -  6.25.  For  example,  K,  run  pair  1  consists  of 
runs  2  and  1  from  Table  D.2.  Umax  run  pair  1  consists  of  runs  33  and  1,  and  y-max 
run  pair  2  consists  of  runs  34  and  2  from  Tables  D.2  and  D.3. 

Table  D.l:  Values  of  Biological  Parameters  Used  in  BI02D  Individual  Sensitivity 
Runs  (see  Table  2.2  for  parameter  definitions). 


RUN 

Mm  a.* 

K. 

Ka 

K, 

Y 

’*6 

Ai, 

F 

Cc 

kc 

1 

6.36 

49.6 

1 

356.8 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

3 

15.36 

49.6 

1 

356.8 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

3 

1.97 

49.6 

1 

356.7 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

4 

6.48 

366 

1 

356.8 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

5 

6.48 

1 

1 

356.8 

.70 

.05 

15.26 

1.0 

3 

715 

,00001 

6 

6.48 

49.6 

1 

1463 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

7 

6.48 

49.6 

1 

23 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

8 

6.48 

49.6 

356.8 

1.02 

,05 

15.26 

1.0 

3 

715 

.00001 

9 

6.48 

49.6 

1 

356.8 

.50 

.05 

15.26 

1.0 

3 

715 

.00001 

10 

6.48 

49.6 

1 

356.8 

.70 

.10 

15.26 

1.0 

3 

715 

.00001 

n 

6.48 

49.6 

1 

356.8 

.70 

.001 

15.26 

1.0 

3 

715 

.00001 

13 

6.48 

49.6 

3 

356.8 

.70 

.05 

15.26 

1.0 

3 

715 

.00001 

13 

6.48 

49.6 

.1 

356.8 

.70 

.05 

15.26 

1.0 

3 

716 

,00001 

14 

6.48 

49.6 

1 

356.8 

.70 

,05 

22.97 

1.0 

3 

715 

.00001 

IS 

6.48 

49.6 

1 

356.8 

.70 

.05 

7.55 

1.0 

3 

715 

.00001 

16 

6.48 

49.6 

1 

356.8 

.70 

.05 

15.26 

1.1 

3 

715 

.00001 

ir 

6.48 

49.6 

1 

356.8 

.70 

.05 

15.26 

.9 

3 

715 

.00001 

i 


Table  D.2:  Values  of  Biological  Parameters  Used  in  BI02D  For  Combined  Sensitivity 
Runs  1-32  (see  Table  2.2  for  parameter  definitions). 


RUN 

K, 

A-„ 

Y 

K, 

F 

Cc 

1 

1.97 

1 

.001 

.5 

.001 

7.55 

356.7 

3 

715 

.00001 

2 

1.97 

266 

.001 

.5 

.001 

7.55 

356.7 

3 

715 

.00001 

3 

1.97 

2 

.5 

.001 

7.55 

356.7 

3 

715 

.00001 

4 

1.97 

266 

2 

.5 

.001 

7.55 

356.7 

3 

715 

.00001 

5 

1.97 

1 

.001 

1.02 

.001 

7.55 

356.7 

1 

3 

715 

.00001 

6 

1.97 

266 

.001 

1.02 

.001 

7.55 

356.7 

1 

3 

715 

.00001 

7 

1.97 

1 

2 

1.02 

.001 

7.55 

356.7 

1 

3 

715 

.00001 

8 

1.97 

266 

2 

1.02 

.001 

7.55 

356.7 

1 

3 

715 

.00001 

9 

1.97 

1 

.001 

.5 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

10 

1.97 

266 

.001 

.5 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

n 

1.97 

1 

2 

.5 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

12 

1.97 

266 

2 

.5 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

13 

1.97 

1 

.001 

1.02 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

14 

1.97 

266 

.001 

1.02 

.1 

7.55 

356.7 

1 

3 

715 

,00001 

15 

1.97 

2 

1.02 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

16 

1.97 

266 

2 

1.02 

.1 

7.55 

356.7 

1 

3 

715 

.00001 

17 

1.97 

.001 

.5 

.001 

22.97 

356.7 

1 

3 

715 

.00001 

18 

1.97 

266 

.001 

.5 

.001 

22.97 

356.7 

1 

3 

715 

.00001 

19 

1.97 

1 

2 

.5 

.001 

22.97 

356.7 

1 

3 

715 

.00001 

20 

1.97 

266 

2 

.5 

.001 

22.97 

356.7 

1 

3 

715 

.00001 

21 

1.97 

1 

.001 

1.02 

.001 

22.97 

356.7 

1 

3 

715 

.OOOOl 

22 

1.97 

266 

.001 

1.02 

.001 

22.97 

356.7 

1 

3 

715 

.00001 

23 

1.97 

1 

2 

1.02 

.001 

22.97 

356.7 

I 

3 

715 

.00001 

24 

1.97 

266 

2 

1.02 

.001 

22.97 

356.7 

J 

3 

715 

.00001 

25 

1.97 

1 

.001 

.5 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

26 

1.97 

266 

.001 

.5 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

27 

1.97 

1 

2 

.5 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

28 

1.97 

266 

2 

.5 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

29 

1.97 

1 

.001 

1.02 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

30 

1.97 

266 

.001 

1.02 

.1 

22.97 

356.7 

1 

3 

715 

.OOOOl 

31 

1.97 

1 

2 

1.02 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

32 

1.97 

266 

2 

1.02 

.1 

22.97 

356.7 

1 

3 

715 

.00001 

168 
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Table  D.4:  Run  Numbers  Comprising  BI02D  Combined  Sensitivity  Run  Pairs  Shown 
in  Figures  6.20  -  6.25.  For  example,  run  pair  1  consists  of  runs  2  and  1,  and  fimaz 
run  pair  1  consists  of  runs  33  and  1. 


Runs  from  Tables  D.2  -  D.3  In  Pair 

PAIR 

K. 

t^maa 

Ko 

Y 

1‘t, 

Kl, 

1 

2-  1 

33-  1 

3-  1 

5-  1 

9-  1 

17-  1 

2 

4-  3 

34-  2 

4-  2 

6-  2 

10-  2 

18-  2 

3 

6-  5 

35-  3 

7-  5 

7-  3 

11-  3 

19-  3 

4 

8-  7 

36-  4 

8-  6 

8-  4 

12-  4 

20-  4 

s 

10-  9 

37-  5 

11-  9 

13-  9 

13-  5 

21-  5 

6 

12-11 

38-  6 

12-10 

14-10 

14-  6 

22-  6 

7 

14-13 

39-  7 

15-13 

15-11 

15-  7 

23-  7 

8 

16-15 

40-  8 

16-14 

16-12 

16-  8 

24-  8 

9 

18-17 

41-  9 

19-17 

21-17 

25-17 

25-  9 

10 

20-19 

42-10 

20-18 

22-18 

26-18 

26-10 

11 

22-21 

43-11 

23-21 

23-19 

27-19 

27-11 

13 

24-23 

44-12 

24-22 

24-20 

28-20 

28-12 

13 

26-35 

45-13 

27-25 

29-25 

29-21 

29-13 

14 

38-27 

4&-14 

28-26 

30-26 

30-22 

30-14 

15 

30-29 

47-15 

31-29 

31-27 

31-23 

31-15 

16 

32-31 

48-16 

32-30 

32-28 

32-24 

32-16 

17 

34-33 

49-17 

35-33 

37-33 

41-33 

49-33 

18 

36-35 

50-18 

36-34 

38-34 

42-34 

50-34 

19 

38-37 

51-19 

39-37 

39-35 

43-35 

51-35 

20 

40-39 

52-20 

40-38 

40-38 

44-36 

52-36 

21 

42-41 

53-21 

43-41 

45-41 

45-37 

53-37 

22 

44-43 

54-22 

44-42 

46-42 

46-38 

54-38 

23 

4&-45 

55-23 

47-45 

47-43 

47-39 

55-39 

24 

48-47 

56-24 

48-46 

48-44 

48-40 

56-40 

25 

50-49 

57-25 

51-49 

33-49 

57-49 

57-41 

26 

52-51 

58-28 

52-50 

54-50 

58-50 

58-42 

37 

54-53 

59-27 

55-53 

55-51 

59-31 

59-43 

28 

56-55 

60-28 

56-54 

56-52 

60-52 

60-44 

29 

58-57 

61-29 

59-57 

61-57 

61-53 

61-45 

30 

60-59 

62-30 

60-58 

62-58 

62-54 

62-46 

31 

62-81 

63-31 

63-81 

63-59 

63-55 

63-47 

32 

64-63 

64-32 

64-62 

64-60 

64-56 

64-48 
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